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Preface 

The  purpose  of  this  study  was  to  develop  a  computer 
code  which  would  validate  the  use  of  the  Boundary  Element 
Method  (BEM)  in  problems  of  fracture  mechanics.  Specifi¬ 
cally,  I  used  the  Displacement  Discontinuity  Method  (devel¬ 
oped  by  Stephen  L.  Crouch)  to  solve  several  fracture 
conf  igur  at  ions . 

Originally,  I  modified  the  code  only  slightly  in  order 
to  solve  problems  involving  isotropic  materials.  Having 
obtained  results  for  those  cases,  I  then  modified  the  code 
even  further  to  include  specially  orthotropic  materials 
(such  as  0°  and  90°  composite  laminates). 

For  introducing  me  to  the  concepts  of  the  BEM,  and 
especially  for  his  patience  and  understanding  during  these 
past  several  months,  I  thank  my  faculty  advisor.  Dr.  Anthony 
Palazotto.  For  coming  up  with  the  displacement  disconti¬ 
nuity  solution  for  am  orthotropic  material,  I  give  credit  to 
divine  intervention.  Perhaps  most  importantly,  for  helping 
to  break  the  tension  which  seemed  to  build  to  new  levels 
each  week,  I  am  grateful  for  the  presence  of  all  the  regular 
participants  in  “Friday  at  the  FlyWright.** 

Ralph  E.  Urch 
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This  investigation  presents  analyses  of  several  frac¬ 
ture  mechanics  problems  via  the  Boundary  Element  Method. 
Specifically,  am  indirect  procedure  known  as  the  Displacement 
Discontinuity  Method  was  used  to  solve  problems  involving 
cracks  in  isotropic  or  specially  orthotropic  materials. 

Infinite  as  well  as  finite  regions  were  considered. 

A  series  of  configurations  were  analyzed  and  compared 
with  either  analytic  solutions  or  results  from  a  finite 
element  model.  Agreement  for  the  inf  ini  te- domain  problems 
was  excellent,  while  solutions  to  the  finite- domain  problems 
ranged  from  good  to  excellent. 

Advantages  and  disadvantages  of  the  Displacement  Dis¬ 
continuity  Method  are  briefly  discussed.  The  main  advantage 
of  the  method  is  the  requirement  to  model  only  the  boundary 
of  the  problem  under  consideration.  The  major  disadvantage 
is  the  time  required  to  solve  the  resulting  fully- populated 


matrix  equation. 

Separate  FORTRAN  codes  are  provided  as  appendices  for 
the  two  material  types  —  isotropic  and  orthotropic.  These 
programs  may  be  utilized  for  either  stress  or  fracture 

analyses.  Program  outputs  include  displacements,  stresses, 

s- 

and  stress  intensity  factors,  as  appropriate.  I  _ 

i 

A 


X 


VALIDATION  OF  THE  BOUNDARY  ELEMENT  METHOD  APPLIED  TO  COMPLEX 
FRACTURE  MECHANICS  CONDITIONS 

I.  Introduction 

In  analyzing  aircraft  structures,  a  mathematical  model 
is  usually  developed  which  involves  an  approximate  represen** 
tation  of  the  physical  component  under  consideration.  The 
most  common  methods  discretize  the  structure  into  a  manage- 
able  number  of  smaller  elements  for  which  a  mathematical 
solution  is  available.  These  methods  may  be  divided  into 
two  groups  Csee  Figure  1.1)  based  upon  whether  the  entire 
domain  or  Just  the  boundary  is  modeled. 

The  most  common  domain  technique  is  the  Finite  Element 
Method  CFEM),  “which  is  a  piecewise  variation  to  minimize 
the  total  potential  energy”  <11:1)  of  the  system.  Another 
example  is  the  Finite  Difference  Method,  in  which  the 
differential  equations  of  the  system  are  approximated  by 
finite  differences  C Newton  forward  difference,  backward 
difference,  central  difference,  etc.).  In  all  cases,  the 
domain  techniques  are  characterized  by  an  approximation  of 
the  system  partial  differential  equations.  The  boundary 
conditions,  however,  are  satisfied  exactly.  (See  Cook  (5) 


and  Desai  <9).) 


Figure  1,1.  Numerical  Met.  bode  of  St.ruct.ural  Analyeie 

On  the  other  hand,  boundary  techniques  approximate  the 
boundary  conditions  of  th«  problem  but.  they  exactly  satisfy 
the  system  partial  differential  equations.  These  methods 
can  be  further  classified  into  direct  or  indirect  methods, 
according  to  whether  the  boundary  parameters  are  solved  for 
directly  or  from  some  other  system  parameters  C2:2>. 

The  most  common  direct  boundary  element  technique  is 
the  Boundary  Integral  Equation  (BIE>,  or  Boundary  Element 
Method  CBEM).  This  involves  combining  the  problem  at  hand 
with  another  <for  which  the  solution  is  known)  in  such  a  way 
that,  when  the  strain  energy  of  the  system  is  integrated 
over  the  volume,  Breen's  formula  may  be  used  to  translate 
the  volume  integral  into  integrals  involving  only  the 
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boundary  of  the  problem.  Applying  the  required  boundary 
conditions  then  yields  the  solution.  (See,  for  example, 

Rizzo  (16),  Snyder  and  Cruse  (21),  Fenner  <10>,  Rizzo  and 
Shippy  C17>,  Cruse  (8),  Okada  et  al  (14),  Shi  and  B ezine 
(19),  and  Shih  and  Palazotto  (20).) 

A  lesser-known  direct  method,  called  the  Fictitious 
Stress  Method  (FSM),  discretizes  the  boundary  into  a  number 
of  straight  line  segments  over  which  the  traction  is  assumed 
constant.  As  will  be  shown  later,  it  is  related  to  the 
indirect  method  chosen  for  this  study  —  the  Displacement 
Discontinuity  Method  (DDM)  —  both  of  which  are  presented  by 
Crouch  and  Starfield  (7). 

The  DDM  is  well-suited  to  problems  in  fracture 
mechanics  because  it  uses  as  its  basic  element  a  line  crack. 
That  is,  each  element  is  considered  to  be  composed  of  two 
surfaces  which  are  coincident  along  a  line  segment.  By  a 
judicious  selection  of  the  definitions  of  each  element, 
problems  involving  cracks  in  a  finite  or  in  an  infinite  body 
may  be  solved. 

Besides  being  able  to  provide  solutions  to  problems 
involving  infinite  domains,  the  boundary  techniques  have  an 
advantage  over  the  domain  techniques  due  to  the  requirement 
to  model  only  the  boundary.  This  reduction  of  order  of  the 
problem  is  especially  useful  considering  that  fewer  elements 
(and  therefore  fewer  degrees  of  freedom)  must  be  used  to 
solve  a  given  problem. 
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There  are  drawbacks,  however.  The  domain  techniques 
are  able  to  take  full  advantage  of  efficient,  banded-matrix 
solution  routines.  The  boundary  techniques,  on  the  other 
hand,  produce  fully- populated  matrices,  which  increases  the 
computational  time  required  for  a  given  number  of  degrees  of 
freedom.  Thus,  there  is  a  trade-off  between  both  techniques 
in  problems  requiring  a  large  number  of  elements. 

Nevertheless,  a  boundary  element  technique  — 
specifically,  the  Displacement  Discontinuity  Method  —  was 
used  in  this  study  and  its  use  in  problems  of  fracture 
mechanics  validated.  First,  a  general  outline  of  the  method 
is  presented,  followed  by  its  use  for  an  isotropic  material. 
Next,  application  of  the  DDM  to  problems  involving  a 
specially  orthotropic  material  is  demonstrated.  For  both 
cases,  stress  analysis  and  fracture  mechanics  applications 
are  reviewed.  Finally,  some  of  the  advantages  and 
disadvantages  of  the  DDM  are  briefly  discussed. 

Included  as  appendices  are  the  two  FORTRAN  programs 
used  in  this  study.  The  first  program  CTVODD  —  Appendix  C> 
was  extracted  directly  from  Crouch  and  Star  field  <7:293- 300> 
for  the  analysis  of  an  isotropic  material.  This  program  was 
modified  only  slightly  to  include  computation  of  stress 
intensity  factors.  The  second  program  CTWODDO  —  Appendix 
D)  is  a  modification  of  the  former  to  incorporate  a 
specially  orthotropic  material. 
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II.  Theoretical  Discussion 


The  DDM  involves  the  concept  of  a  displacement  discon¬ 
tinuity  along  a  line  crack,  which  may  be  imagined  "as  a  line 
crack  whose  opposing  surfaces  have  been  displaced  relative 
to  one  another."  <7:79)  The  DDM  uses  the  special  case  where 
the  relative  displacements  are  constant  along  the  entire 
crack.  In  reality,  the  relative  displacements  vary  over  the 
entire  length  of  the  crack;  however,  if  the  line  crack  is 
discretized  into  a  sufficient  number  of  smaller  line  cracks 
< elements),  the  displacements  will  be  approximately  constant 
over  the  length  of  each  element.  The  solution  obtained  by 
superposition  of  the  displacement  fields  of  all  the  elements 
will  then  represent  the  true  solution  with  sufficient  accu¬ 
racy.  This  is  the  basis  of  the  DDM. 

If  a  series  of  line  cracks  in  an  infinite  plane  are 
Joined  end-to-end  to  form  a  closed  contour,  two  regions  are 
formed.  The  interior  region  may  be  used  to  obtain  the  solu¬ 
tion  to  a  problem  involving  finite  dimensionality,  while  the 
exterior  region  may  be  used  to  represent  a  cavity  in  an 
infinite  body. 

As  will  be  seen  later,  the  DDM  involves  the  solution  of 
a  system  of  equations  which  are  determined  by  specifying 
either  <i)  the  tractions  on  each  element,  <ii>  the  displace- 
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merits  of  one  of  the  surfaces  of  each  element,  or  Ciii>  a 
combination  of  the  two. 

A.  The  Method 

The  Displacement  Discontinuity  Method  is  characterized 
by  use  of  the  solution  to  a  line  crack  in  an  infinite  plate 
where  the  two  surfaces  of  the  crack  have  a  constant  displacement 
with  respect  to  one  another  CFigure  2.1  >.  The  displacement 
field  for  an  arbitrarily  oriented  crack  <element>  is 
obtained  from  this  solution  via  coordinate  transformation. 

For  two  or  more  elements,  linear  superposition  is  used 
by  alternately  considering  only  one  element  to  have  a  dis¬ 
placement  discontinuity  present  <i.e.  —  the  displacements 
are  considered  continuous  at  the  other  element  locations). 

The  sum  of  each  of  these  displacements  then  yields  the  dis¬ 
placement  due  to  the  presence  of  all  the  elements,  assuming 
of  course  that  the  material  is  linearly  elastic. 

By  applying  appropriate  boundary  conditions,  a  linear 
system  of  equations  is  formed.  Solution  of  this  system  for 
the  relative  displacements  at  each  element  then  allows 
determination  of  the  displacements  and  stresses  at  each 
element.  The  stresses  and  displacements  at  other  points  in 
the  body  can  be  obtained  based  on  these  same  relative  sur¬ 
face  displacements. 

The  procedure  may  be  divided  into  the  following  11 
steps: 
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Requirement.: 

Displacements  are  continuous  everywhere 
except  across  the  crack. 

Normal:  u  -  u+  -  D  *  constant 
n  n  n 

Shear:  u  -  u  ■  D  ■  constant 

9  9  9 


Figure  2.1.  Line  Crack  in  an  Infinite  Plate 


<1>  Use  the  solution  to  the  problem  of  a  constant  displace 
ment  discontinuity  along  a  line  segment  <crack>  in  an  infi¬ 
nite  plate,  which  takes  the  following  form: 


ds 

/,!«.?>  Ds 

Ds 


+  Dn 

♦  V*'y>  °n 

*  Dn 
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a_  > 
yy 


C2.4) 


'  ts®*’  D.  ♦  Dn 

Txy*  /5x<Sl5i>  °s  *  *„<**>  Dn  <2'55 

wh«r«  x,  y  Cor  s,  n)  is  ths  local  coordinates  system  Csss 

Figure  2.2),  the  functions  /  _  and  f  _  depend  only  upon  the 

local  coordinates  and  the  length  of  the  crack,  and  D  and  D 

s  n 

are  the  shear  and  normal  displacement  discontinuities, 
respectively.  That  is,  is  a  (displacement)  discontinuity 
in  the  x-direction  while  is  a  discontinuity  in  the 
y- direction  such  that 


C2.6) 

C2.7) 


and 


<«•  »» 


refer  to  the  positive  and 


where  the  superscripts  *-.v. 
negative  surfaces,  respectively,  of  the  crack. 


C2>  Use  coordinate  transformation  to  relate  the  local  coor¬ 
dinates  Cx,  y>  to  the  global  coordinates  Cx,  y)  for  an  arbi¬ 
trary  rotation  CJ/£3>  and  translation  CJx,  Jy>  as  depicted  in 
Figure  2.2: 


x  ■  Cx  -  Jx>  cos  }0  *  Cy  -  sin  *0  <2.8> 

y  «  — Cx  -  Jx>  sin  30  +  Cy  -  Jy>  cos  30  C2.9> 


u  * 

X 

u 

X 

cos  J0 

-  u_  sin 

y 

C2.10> 

u  ■ 

y 

u_ 

X 

sin 

♦  u_  cos  JS 
y 

C2.11> 

2 

a  «■  cr  _  cos 

XX  XX 

*0 

-  r 

xy 

sin  2J£  cy _ 

yy 

sin2  30 

C2.12> 

a  m  a _ sin2 

yy  xx 

1 0 

+  T _ 

xy 

sin  2*0  +  o _ 

yy 

cos2  30 

C2.13> 

*•  (<? _ -  a _ ) 

xx  yy 

sin  3 0  cos  30  *  r _ Ceos2  30  -  sin2 

xy 

j0> 

C2.14> 


C3>  Consider  a  series  of  N  line  cracks  Celements>  joined 
end-to-end  CFigure  2.3>.  Determine  the  effect  on  the  ith 
element  of  displacement  discontinuities  at  the  jth  element. 
Because  the  displacements  and  stresses  prescribed  at  each 
element  are  specified  in  the  local  coordinate  directions  of 
the  element,  transform  the  above  equations  to  the  tth 
element’s  local  system  Cx’,  y’>: 


'’u.  ■  u  cos  iB  u  sin  i0  C2.15> 

x'  x  y 

lu.  ■  -u  sin  i0  u  cos  ^ 0  C2.16> 

y'  x  y 

vO'_  _  ■  a  com2  l0  *  t  sin  2i0  ♦  a  sin2  *0  C2.17> 

k'x'  xx  xy  yy 
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X 


Figure  2.3.  S«rlM  of  N  Lino  Crack*  Joined  End- to- End 

vo' _ ■  o  sin2  l0  -  t  sin  2^0  +  a  co*2  V  <2.18> 

y'y'  xx  xy  yy 

Vr _ ■  -(.o'  -or  >  *in  coop  *  r  <co*2  ^  -  sin2  *'fl) 

x'y'  xx  yy  xy 

C2.19) 

<4>  Not.*  that  th*  displacanwni*  vu_  and  u_  and  tha 

X'  y' 

stresses  lo_  and  vt_  _  are  equivalent  to  the  local  shear 

y'9  x'y' 

and  normal  displacements  and  stresses  at  the  ith  element: 
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(2.20) 


L 

O' 

s 


l 

cr 

n 


V 

T 

x'y' 

l 

a _ 

y'y' 


(2.21) 

(2.22) 

(2.23) 


(5)  Use  the  above  four  steps  to  express  the  effect  of  dis¬ 
placement  discontinuities  JD  and  JD  at  the  jth  element  on 

s  n 

the  displacements  and  stresses  at  the  center  of  the  ith 
element: 

lu  ■  u  cos  V  +  u  sin  (2.24) 

S  x  y 

lu  -  -  u  Sin  ^0  u  cos  (2.25) 

n  x  y 

o'  ■  “(o'  “O'  )  sin  cos£  +  r  (cos2  -  sin2  l0) 

S  xx  yy  xy 

(2.26) 

"o'  -  o-  sin2  -  r  sin  2V0  +  o-  cos2  (2.27) 

n  xx  xy  yy 

where  subscripts  “x‘*  and  "y"  refer  to  the  global  coordinate 
directions,  and 


u  ■  u_  cos  20  -  u_  sin  J0 
xx  y 

■  /t _  (S , y) J D^co*  V  ♦  /i_Cx,y)  JDncos  V 
-  /2S<*»y>JI>ssin  V  “  /a_(x,y)JDnsin 

■  j/lS<S,y>  cos  }0  -  /2_(x, y>  sin  joJjDs 

♦  £/t-(x,p)  cos  j0  -  /2-<S,p)  sin  VjjDn  (2.28) 
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u  ■  u_  sin  3  0  +  u  cos  J0 
y  x  y 

■  /.-<x,y>JD  sin  30  *  /  _<x,y>JD  sin  30 


iy 


♦  f2  _<x,y>  JDflcos  ■*  /2_<x,y>  ^D^cos  30 

[/  _<x,y>  sin  3  0  f  _<x,y>  cos  J0| '’d 

1  ^  2  x  I  $ 

♦  ^/^-CXjp)  Sln  +  -^2y**’^  COS 


o'.  _  cos2  30  -  t _ sin  230  +  o _ sin2  30 

xx  xy  yy 

/  _<x,y>^D  cos2  3  0  +  f  _<x,y>^D  cos2  ^ 

3  X  &  3  y  R 

-  /  _(x,y>^D  sin  2J£  -  /  _Cx,y>*D  sin  2^ 

3x  s  ■'sy  '  n 

+  /  _Cx,y>,’D  sin2  3  0  +  f  _<x,y>^D  sin2  30 

4X  7  S  J 4y  z  n 

[/  _<x,y>cos2  3 0  -  f  _<x,y>sin  2J£ 

3  x  3  x 


+  /  _<x,y>sin 

4x  7 


J«] ’D.  +  [',*<*■ 


y>cos2 


yy 


-  /5_<x,y)sin  230  +  /^_<x,y>sin2 
*  a _ _  sin2  t _ sin  230  +  a _ cos2  30 

xx  xy  yy 

"  /3i«'y>JVin2  *0  *  /#_<x,y>jDnsin2  V 

+  /  _<x,y> JD  sin  23 0  +  /  _<x,y>JD  sin  2J0 
5x  s  ■'sy  '  n 

+  /^_<x,y>  ^Dgcos2  ♦  /^_Cx ,y>  ■’d^cos2  j0 


y>sin2  30  f  _ Cx , y>sin  2J/9 


Sx 


/  _<x,y>cos2  ^"1  +  [/  _<x,y>sin2  30 

4X  J  ^  L 


♦  /__<* ,y>»in  aj£  /^_<x,y>oos2  Jpj  Jl 


Sy 


n 


<2.29> 


<2.30) 


<a.ai> 
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T*v  ■  <a55  “  °'cc>  sin  cos  +  T--  <cos2  J£-  sin2  J0) 

xy  xx  yy  xy 

•  |/3x<s'y>jDs  +  /39<s>y>jDn  •  ^s<s'y)Jl>8 

-  /  _<x,y)JD  sin  30  cos  30  +  |/  _<x,y)^D 
4y  n_  [_  Sx  s 

+  /  _<x,y)^D  (cos2  -  sin2  3 0) 

5y  nj  ^  ^ 

-  |[/3S<S,y>  -  /45<x,y>J  Sin  V  cos  J0 
+  /s_<x,y)<cos2  -  sin2  J£)^De 

♦  "  /4,<-x,y>]  sin  J0  cos 

+  /5_<x,y)(cos2  V  -  sin2  J0>|JDn  C2.32> 

By  properly  defining  a  new  set  of  functions  /  and 
/  ^,  the  global  displacements  and  stresses  can  be  written  in 


the  form 

u  -  /  <x,y,j£)jD  ♦  /  <x,y,J0>JD  C2.33> 

X  •'lx  s  ly  7  n 

u  -  /  <x,y,j£)jD  +  /  <x,y,^0)^D  <2.34) 

y  2x  S  2y  n 

cf  rn  /  <x,y,J0)JD  ♦  /  <x,y,J0y*D  C2.35> 

XX  J  3x  S  ^  3y  7  K  n 

O'  -  /  <x,y,j0)JD  +  /  <x,y,J0)JD  C2.36> 

yy  S  4y  n 

r  «  /  <x,y,J/9)jD  ♦  /  <x,y,j0)jD  <2.37) 

xy  "'Sx  7  K  S  3y  7  K  n 


while,  by  using  Eqs  <2.24)  through  <2.27),  the  local  dis 
placements  and  stresses  take  the  following  form: 


Vu  -  VJB  <x,y,J0)JD  +  vjB  Cx,?,^)^  <2.38) 

s  as  s  sn  n 

vu  -  VJB  <x,y,J^)jD  +  vjB  <x,y,V>JD  <2.39) 

n  ns  s  nn  K  n 

X*m  -  ijAss<R,y,V>jDs  ♦  ijAsn<x,y,V>jDn  <2.40) 

lo'n  -  ijAr-<R,y,i0)jD-  ♦  ^A^x,?,'*)^  <2.41) 
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where  lJB  ,  ',JB  ,  ...»  lJA  are  referred  to  as  influence 
ss  sn  nn  J 

coefficients. 

C6>  Compute  the  total  effect  at  the  ith  element  as  the  sum 
of  the  effects  of  all  N  elements: 


re  n 

\  lJB  jD  +  \  ljB  JD 

/  ss  s  /  sn  n 

j  =  i  j=  i 

N  N 

S  VJB  JD  ♦  \  VJB  jD 

/  ns  s  /  nn  n 

j  =  i  j  =  i 

N  N 

\  tJA  jD  +  \  VJA  JD 

_ _  ss  s  /  sn  n 

j  =  i  j=i 

H  N 

V  VJA  >  ♦  V"  ijA  JD 

/  ns  s  /  nn  n 


<2.42> 


C2.43> 


<2.44> 


<2.45> 


C7>  Determine  the  effect  of  the  ith  element  on  itself 
(element  self-effects >  by  computing  the  influence  coeffi¬ 
cients  for  the  case  that  i  ■  j  and  x  *  y  «  0.  The  result  is 
that  all  of  the  stress  influence  coefficients  c'-lA  's>  are 
single-valued,  but  that  some  of  the  displacement  influence 
coefficiente  <UB  '»>  are  dual-valued.  Namely, 


-  UB  -  +-  for  y  “  0  <2.46> 

ss  nn  a  ± 

The  discontinuous  nature  of  VVB  and  '’>>B  does  not 

ss  nn 

cause  complications,  however,  if  a  direction  of  travel  is 
chosen  consistently  when  defining  the  boundary  of  a  problem. 
Referring  to  Figure  2.4,  Crouch  and  Starfield  C7:64-63> 
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Direction  of  trave*  shown  is  for  an  exterior 
problsm  (cavity  in  an  infinite  plane). 


O 

Figure  2.4. 


Defining  the  Boundary  of  a  Problem 


state  that 

the  values  of  coefficients  £  B8fl]  and  £  depend 

upon  the  way  in  which  curve  C'  is  approached.... 

Curves  CJ.  and  can  be  regarded  as  the  boundaries  of 
the  interior  and  exterior-  rofions,  respectively,  of 
curve  C'....  tin  we  traverse  curve  C'  in  the 
clockwise  sense  ...  then  the  negative  side  of  the  curve 
defines  the  boundary  of  the  interior  region.  We  can 
therefore  avoid  having  to  use  different  values  for  the 
coefficients  t  1  and  £  *  B  J  if  we  agree  slwstys  to 

traverse  the  boundary  in  such  a  way  that  the  outward 
normal  points  away  from  the  region  of  Interest. 
Accordingly,  we  adopt  the  following  convention:  The 
boundary  of  a  finito  body  is  troversecf  in  the  clockwiso 
sense,  whor&as  the  boundary  of  a  cavity  is  traversecf  in 
the  countorclockwiso  sense.  This  convention  leads  to 
some  simplification  in  comp< 
means  that  coefficients  £  i 
to  *  1/2  for  both  types  of 


iter  programming  because  it 
1  and  rvB  1  are  equal 
^Toblems.  ^ 


2- 


<8>  Choose  the  equations  appropriate  for  the  prescribed 
boundary  conditions  of  the  problem.  The  boundary  conditions 
at  each  element  can  be  grouped  into  four  categories: 


<i  ) 

i. 

cr 

s 

and 

V 

cr 

n 

are 

prescribed 

<ii> 

l 

u 

8 

and 

i 

u 

n 

are 

prescribed 

<ii  i  > 

l 

u 

s 

and 

i 

a 

n 

are 

prescribed 

<iv> 

l 

Cr 

s 

and 

i. 

U 

n 

are 

prescribed. 

Thus,  a  system  of  2N  equations  is  formed: 


N 


N 


lJc  JD  i  -  1,2,...N 

sn  n 


lb  -  S  lJc  JD 

s  _ _  ss  s  /  . 

j= i  j  =  i 

N  N 

vb  -  \  VJc  JD  +  VJc  JD  i  -  1,2, ...N 

n  /  ns  s  ^ _ _  nn  n 


j  =  i 


j*  i 


where 


<2.47> 


C2.48> 


v  W  —  L 

y  or  u 

b  - 

V  1 

cr  or 

u 

& 

s 

s 

n 

n 

n 

VJc 

-  VJA 

or 

VJB 

VJc 

< 

*n 

■ 

or 

■0» 

ss 

ss 

ss 

ns 

ns 

ns 

VJC 

-  V  JA 

or 

VJB 

vjc 

< 

n 

a 

or 

VJB 

sn 

sn 

sn 

nn 

nn 

nn 

<2.49> 


as  appropriate. 

<9>  Solve  the  system  of  2N  equations  for  the  2N  unknowns  CN 

unknown  JD  N  unknown  JD  -  2N  unknown*)  using  direct 
s  n 

Gaussian  elimination.  Because  of  the  diagonal  dominance  of 
the  matrix,  no  pivoting  is  required  to  ensure  an  accurate 
result. 


2-12 


<10>  Compute  the  displacements  and  stresses  at  each  element 
using 


N 


N 


U 


U 


n 


n 


S  VJB  JD  ♦  \  lJB  3D 

_  ss  s  /  sn  n 

j  =  i  j=  i 

N  N 

\  VJB  3D  +  \  lJB  JD 

Z_  08  s  Z_ 

j  =  i 

_  N 

\  tJA  JD  +  \  lJA  JD 

_ _  ss  s  /  sn  n 

J  =  »  j=  i 

N  N 

VJA  JD  +  \  '  VJA  JD 

_ _  ns  s  ^ _ _  nn  n 


ns  s  ^ _ _  nn  n 

j  =  i  j  =  i 

N  N 


J  =  1 


J=  1 


l  ♦  V  - 

u  ■  u  -  D 

8  S  S 


i.  +  v  -  i_ 

u  ■  u  -  D 
n  n  n 


<11  >  Compute  displacements  and  stresses  at  other  points 
the  body  using 


N 


U 


U  a 


‘Dn 


yy 


j  =  i 
N 

E 

j  =  i 

N 

E 

j- 1 
N 

E 

j=  i 


/2x<S,y,J0>  JDS 


/3x<x, y,J/9>  JDt 


E 

j= i 
N 

7""  /,„<x,y,^>  JDr 


2y 


J=  l 
N 


E  V8*^ 

J  =  l 
N 

E 


♦  )  /4y<S,y,J0>  JDn 

j=  i 


<2.50> 

<2.51> 

<2.52> 

<2.53> 

<2.54> 

<2.55> 

in 

C2.56> 

<2.57> 

<2.58) 

<2.59) 


2*13 


<2.60> 


N  N 

T,y  ’  XI  \  *  XI  V5'5^5  'Dn 

J=  i  J  =  1 

where 

5  -  Cx  -  Jx>  cos  J0  Cy  -  V>  sin  *0  C2.61> 

y  -  -<x  -  Jx>  sin  J0  +  Cy  -  Jy>  cos  (2.62) 


and  (Jx,Jy)  are  the  midpoint,  coordinat.es  or  the  jth  element.. 


B.  Computer  Implementation 

The  DDM  lends  itself  to  being  programmed  into  a  comput¬ 
er  with  ease  due  to  the  modular  nature  of  the  procedure. 
Appendix  C  provides  a  listing  of  the  FORTRAN  code  for  pro¬ 
gram  TWODD,  along  with  an  explanation  of  the  input  deck. 

First,  the  material  properties  and  symmetry  conditions  are 
defined.  Next,  as  the  boundary  conditions  of  the  problem 
are  input,  the  system  of  equations  is  set  up.  Then,  the 
system  is  solved  using  direct  Gaussian  elimination.  The 
stresses  and  displacements  at  each  boundary  element  are  then 
computed,  followed  by  displacements  and  stresses  at  any 
other  points  of  interest.  Other  features  Cinclusion  of  a 
title,  initializing  counters  for  loops,  etc.)  are  added 
where  appropriate. 

For  an  isotropic  material,  the  only  material  properties 
required  for  input  are  Young's  modulus  CE)  and  Poisson's 
ratio  <v).  The  shear  modulus  CG)  is  computed  from  these  by 


2-14 


G 


1  E 

2  1  +  u 


<2.63> 


For  an  orthotropic  material  E  ,  E  >  G  ,  and  v  must  be 

1  2  12  12 

input. 

If  symmetry  of  geometry  and  loading  exist  in  the 
problem,  they  may  be  incorporated  to  reduce  the  number  of 
unknowns  in  the  matrix  equation.  For  example,  suppose  x  ■  0 
is  a  line  of  symmetry  (Figure  2.5>.  Then  the  conditions  at 
element  2  are  determined  by  those  at  element  1: 


u 

m 

- 

u 

D 

m 

- 

s 

s 

s 

s 

u 

m 

+ 

1 

u 

2d 

m 

4- 

n 

n 

n 

n 

i 

2 

i 

a 

m 

+ 

a 

Or 

m 

Cr 

s 

s 

n 

n 

<2.64> 


Figure  2.9.  Line  of  Symmetry 
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In  place  of  the  original  system  of  equations. 


u 

m 

n 

c 

+ 

12  2_ 

C  D 

11 

C  D 

+ 

12  2_ 

C  D 

(2.65) 

s 

ss  s 

SS  s 

sn  n 

sn  n 

u 

m 

11 

c 

12  2_ 

C  D 

•4- 

11  if* 

c  D 

+ 

12  2 
c  D 

(2.66) 

n 

ns 

s 

ns  s 

nn  n 

nn  n 

u 

m 

21 

c 

+ 

22  2_. 
c  D 

+ 

21 

c  D 

+ 

22  2- 
c  D 

(2.67) 

s 

ss 

s 

ss  s 

sn  n 

sn  n 

u 

m 

21 

c 

+ 

22  2_ 
c  D 

c  D 

22  2_ 
c  D 

(2.68) 

n 

ns  s 

ns  s 

nn  n 

nn  n 

a  new  system  of  equations  may  be  formed: 


u  m 

u 

c 

+ 

12 

C 

2d 

+ 

ii 

C 

+ 

12 

C 

2d 

(2.69) 

s 

ss  s 

SS 

s 

sn 

n 

sn 

n 

u  ■ 

11 

c 

+ 

12 

c 

2D 

n 

c 

lD 

12 

c 

2D 

(2.70) 

n 

ns 

s 

ns 

s 

nn 

n 

nn 

n 

2D  -  -  *D  (2.71) 

s  s 

2D  ■  +  XD  C2.72) 

n  n 


By  substituting  the  last  two  equations  into  the  first  two,-.1 
these  may  be  rewritten  as 


u  * 

,ii 
(  c 

12  Or. 

c  )  D 

+ 

,ii 
(  c 

12 

c  )  D 

s 

ss 

ss  s 

sn 

sn  i 

u  m 

,11 
(  c 

12  Or. 

c  )  D 

,n 
(  c 

12  Or. 

c  )  D 

n 

ns 

ns  s 

nn 

nn 

-  ‘D 


D  - 

n 


<2.73> 

(2.74) 

(2.75) 

(2.76) 


Now,  instead  of  solving  four  equations  in  four  unknowns 

<*D  ,  *D  ,  2D  ,  *D  ,),  only  the  first  two  equations  in  two 
s  n  s  n 

unknowns  must  be  solved.  The  conditions  at  element  2  are 
then  computed  from  the  last  two  equations.  Thus,  applica~ 
tion  of  symmetry  cuts  in  half  the  size  of  the  matrix  equa¬ 
tion.  If  two  lines  of  symmetry  exist,  then  the  size  of  the 
matrix  equation  is  reduced  by  a  factor  of  four. 
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In  order  to  solve  problems  involving  infinite  domains 
with  a  remotely  applied  stress,  changes  in  displacement  and 
stress  are  considered.  This  is  done  by  adjusting  the  bound¬ 
ary  conditions  of  each  element  to  include  the  negatives  of 
the  tractions  that  would  exist  if  the  elements  were  not 
present.  The  solution  is  then  obtained  by  summation  of  the 
stress  and  displacement  changes,  assuming  zero  initial  dis¬ 
placement  throughout  the  plane. 

After  the  system  of  linear  algebraic  equations  is  set 
up,  it  may  be  solved  using  direct  Gaussian  elimination  with 
no  pivoting.  This  is  because  the  logarithmic  terms  in  the 
influence  coefficients  make  the  matrix  diagonally  dominant: 
a  nearby  element  has  more  effect  than  a  faraway  element. 
Since  the  diagonal  terms  account  for  element  self-effects, 
these  coefficients  are  larger  in  magnitude  than  the  other 
coefficients  in  the  corresponding  row/column. 

Once  the  discontinuities  at  each  element  are  deter¬ 
mined,  the  stresses  and  displacements  are  computed  readily 
as  the  sum  of  the  influences  from  each  element.  In  order  to 
reduce  the  size  of  the  program,  only  one  subroutine  is  used 
for  all  influence  coefficient  calculations.  The  coeffi¬ 
cients  computed  are  the  functions  /lx>  /ly»  •••*  /5y  in  step 
11  above  CThe  Method).  Coordinate  transformation  of  these 
coefficients  allows  calculating  the  VJA  ’s  and  lJB  'a,  while 
no  transformation  is  required  for  computing  conditions  at 
non- boundary  points. 
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C.  Stress  Intensity  Factors 

When  dealing  with  problems  involving  cracks  in  struc* 
tures,  a  very  useful  parameter  is  the  stress  intensity  fac¬ 
tor,  denoted  by  K.  For  two-dimensional,  elastic  fracture 
mechanics  problems,  this  measure  of  the  intensity  of  the 
stress  singularity  at  the  crack  tip  may  be  decomposed  into 
two  modes  (Figure  2.6):  Mode  I  is  the  “opening"  mode  and 
Mode  II  is  the  “shearing"  mode. 

If  an  infinite  plate  contains  a  single  crack  (Figure 
2.7),  the  stress  intensity  factors  in  a  region  very  close  to 
the  crack  tip  Cthe  near- field  solution )  are  given  by  <3:79, 

81) 

Mode  I: 

K 

o  rn  - 1 cos  0/2  [l  +  sin  B/2  sin  30/2  1  <2.77> 

y  Y-2ZF-  l  J 

Mode  II: 

K 

T  .  - 11—  cos  0/2  |l  -  sin  0/2  sin  30/2  1  <2.78> 

*y  V~2^  J 

For  a  problem  in  which  the  stress  field  in  the  vicinity  of 
the  crack  tip  is  known,  the  stress  intensity  factors  may  be 
determined  by  a  limiting  procedure  (3:343): 

Ki  m  r-»small  ^nr  Tcos  fl  ♦  sin  0/2  sin  30/2 


(2.79) 


MOOC  I: 


O 
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a 

y 


t  r 


Figure  2.7.  Coordinates  for  the  Near-Field  Solution 
Kit  "  r -►small  TxyV  2nr  [cos  0/2  f1  +  ®in  0/2  9in  30/2  )] 

i 

(2.80> 

where  r  is  the  distance  from  the  crack  tip. 

Due  to  numerical  inaccuracies  associated  with  digital 
computers  when  dealing  with  singularities,  the  stress  inten¬ 
sity  factors  in  this  study  were  obtained  by  applying  a  least 
squares  linear  fit  0:264)  to  a  plot  of  anfenr  vs.  r  Cor 
tV2 nr  vs.  r  for  Mode  II).  Here  the  stresses  are  calculated 
directly  in  front  Of  the  crack  C©  w  0).  Figure  2.8  depict® 
a  typical  curve  generated  using  the  DDM,  along  with  the 
associated  linear  fit.  Note  the  Inaccuracy  of  the  solution 
near  the  crack  tip  Cr  ■  0).  Far  from  the  crack  tip,  the 
near- field  solution  is  no  longer  valid.  Therefore,  the 
solution  tends  away  from  a  straight  line  there. 
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TYPICAL  CURVE 


Linear  Fit  between:  0.20  and  0.60 
Intercept  Value:  1.8783 


The  suitability  of  the  straight-line  fit  may  be  deter¬ 
mined  by  examining  the  index  of  correlation  <l:264-287>, 
defined  for  an  nth-order  polynomial  fit  to  m  data  points  by 


/.v 


_  r 


2  1 1/2 


m  o 


)' 


in  which 


o-  -  Y 
Y  var 


-  E- 

k  =  1  t= 1 

m 

■=Ev 


/Li  > 
Y 


J  =  1 


M  -  Y  =  —  y 
Y  av9  in  _ ,  J 

j  =  i 


The  best-fit  polynomial  is 


-  Z  ^ 


V  =  1 


<2.81) 


<2.82) 


<2.83) 


<2.84  > 


C2.85) 


with  the  coefficients  o  being  determined  by  solving  the 
linear  system  of  equations 

b  c  ■  v  i  ■  l,2,...,n;  j  ■  l,2,...,n  <2.86) 

ij  v  j  ~ 

where 


b  *  s 

W  v*j-l 


^IE-c‘ 

k  s  i 


l,2,...,2n-l 


C2.87> 

<2.88> 
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m 


k=  1 


i  •  l,2,...,n 


<2.89> 


and  Cx^,  yfc)  are  the  coordinates  of  the  Ath  data  point. 

For  a  linear  fit,  the  definition  of  the  index  of  corre¬ 
lation  may  be  expressed  as  <12) 


Ev.  -*EnE*. 


k  =  1 


k  =  1  k  =  1 


/.V 


E  ■  s  E  ^  E  <  -  s  E 

k  =  1  k  =  1  k  =  1  k  =  1 


1X2 


<2.90  > 


The  value  of  ryY  niay  range  from  zero  to  unity.  The 
closer  r  is  to  unity,  the  better  the  fit.  A  perfect  fit 

j  ^ 

is  represented  by  a  value  of  exactly  unity,  while  an 
extremely  poor  fit  produces  an  index  close  to  zero. 


D.  Modeling  Considerations 

Because  the  relative  displacements  of  the  surfaces  of  a 
real  crack  are  not  constant  <as  is  assumed  for  the  DDM  ele¬ 
ment),  a  method  of  modeling  a  crack  was  chosen  which  would 
more  closely  resemble  the  true  situation  CFigure  2.9).  Each 
half-crack  (containing  only  one  crack  tip)  was  discretized 
into  at  least  three  DDM  elements.  For  the  coarsest  case, 
the  element  at  the  tip  was  one-sixth  the  length  of  the  half¬ 
crack,  while  the  element  furthest  from  the  tip  was  one- half 
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Physical  Model  of  One  Half-Crack 

crack  tip 


1 

U - a - - 


3 


DDM  Models  of  One  Half- Crack 


Half -Crack  Elements: 


crack  tip 


6  Half-Crack  Elements: 


192  Half-Crack  Elements: 


64  at 


192 


-64  at 


*~64  **•  128 

< -  »  » 

n 

mr  .■  ,i 

l!!'“  •  •  •"•  III 

ii . n 

L  *  r| 

t  a 

L*  j 

^  2 

K  3  * 

^6^ 

2.9.  Modeling  a 

Half- Crack 

Using 

384 


crack  tip 
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the  length  of  the  half-crack.  This  left  the  middle  element 
with  a  length  of  one- third  the  total  length. 

Finer  element  meshes  were  created  by  dividing  each 
element  into  equal  numbers  of  smaller  elements.  Thus,  a 
result  stating  that  “192  half-crack  elements  were  used** 
means  that  64  equally- sized  elements  were  used  in  the  first 
one-sixth  of  the  half-crack,  64  equally-sized  elements  were 
used  in  the  next  one- third,  etc.  This  allowed  the  crack  to 
more  easily  take  on  the  elliptic  shape  it  would  assume  under 
a  tension  loading. 
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III.  Isotropic  Development 


A.  Theoretical  Development 

The  basic  solution  required  for  the  Displacement  Dis~ 
continuity  Method  for  an  isotropic  material  subjected  to 
plane  strain  is  given  by  Crouch  (6>.  Using  the  notation 
from  Chapter  II,  the  displacements  due  to  a  constant  dis~ 
placement  discontinuity  along  a  line  crack  defined  on 
|  x  |  <  a  (Figure  3.1>  is  presented  as 

u_  »  D  [  2(1  -  v>/_-  yf _ 3  +  D  £  -  (1  -  2v>/_-  y/„] 

x  s  ,y  ,xx  n  ,x  ,xy 

(3.1> 

* 


D  ■  u”  -  u+  ■  constant  along  crack 

m  x  s 

D  ■  u"  -  u+  ■  constant  along  crack 

n  y  y 


Figure  3.1.  Displacement  Discontinuity  along  a  Line  Crack 


u_  =  D  I  -  Cl  -  2u~>  f  yf J  +  D  (  2C1  -  *>>/_-  yf  1 

y  S  »X  »xy  n  ,y  '■',55 


<3.2> 


where 
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x  -  a 
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-  arctan 

1X2 


— Y — 1 

x  +  aj 


+  Cx  +  a>  Injcx  +  a>2  +  y 2J 1/2 


—  + 

I  m  U_  -  U_ 

S  x  x 
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on 

1*1 

<  a 

—  ♦ 

1  “  u_  -  u_ 

n  y  y 

■  constant 

on 

1*1 

<  a 

D 
D 

and  the  derivatives  of  /Cx,y>  are 


<3.3> 

<3.4> 

<3.5> 
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p> 
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p> 

+1 
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lxa 


arctan  s--~  “  arctan  ^  ~ 


x  -  a 


Cx  -  a)2  +  y2  Cx  4  a>2  ♦  y2  J 


C3.6> 


C3.7> 


C3.8> 


■^,x5  ^ •yy  4rtCl  -  p> 


x  -  a 


x  ♦  a 


|_  Cx  -  a>2  ♦  y2  Cx  •»•  a>2  ■*•  y2  J 


C3.9> 


The  stresses  for  the  isotropic  material  are  determined 
from  the  stress-displacement  relationships: 


20 

<y  m  ■ 

55  1  -  2u 


Cl  -  p)  u_  ♦  p 


a _ 


2G 


_  i  I  v  u_  _  ♦  <1  -  P> 

yy  1  -  2w  I  x,ti 


u.  1 

y»y  J 

u-l 

y»y  J 


C3.10> 


C3.11> 
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T8y  *  °  (  “5.9  *  “y.5  )  <3  ,2> 

where  G  "  2 OTT)  i®  the  material  sheer  modulus. 

Applying  coordinate  tr ansf ormation,  the  displacemonts 

i 

in  the  global  Cx,  y>  coordinate  system  are  <7:91> 


u  *  D«  "  <*  "  2 v>  sin  0  /_  +  2<1  -  i->>  cos  0  / 
x  ®L  -x  -v 


+  D 


+  9(Sln  e  '.5y  -  008  *  >55)] 

-  <1  -  2v>  cos  0  /_  -  2<1  -  v>  sin  0  /_ 
*L  «x  <9 

-  y(co»  0  /„  ♦  »ln  0  /„]] 

u  ■  D  +  (1  -  2v»>  cos  0  /_  +  2C1  -  sin  0  /_ 

y  ®L  <x  ,y 

-  ?(=<».  0  /55  »  -in  0  /__]] 

♦  Dn[  "  <*  -  2v>  sin  0  /_  +  2<1  -  i>>  cos  0  /_ 

-  y(«n  0  -  00-  e  /„]] 


<3.13> 


C3.14> 


From  these,  the  stresses  in  the  global  coordinate  sys¬ 
tem  are  determined  to  be  <7:92> 

or  ■  20  D  I  2  cos2 0  / _ ♦  sin  20  f 

XX  8  ■'.xy  ,xx 

+  y  [cos  20  / _ -  sin  20  / _ 11 

l.  >xyy  .yyyJJ 

♦  20  D  [-/_.♦  yfsin  20  / _ +  cos  20  / _ H 

n  L  l  J  >xyy  .yyyJJ 


<3.13> 
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®  »xy  ,xx 

-  y  I  cos  20  / _ -  sin  20  / _ H 

l  *xyy  *yyyj  J 


+  2G  D 


n  [  -  '.99  * 


sin  20  / _ ♦  cos  20  / _ 11 

»xyy  .yyyJJ 


<3.16> 
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♦  yf< 


y  |sin  20 


/ _ +  cos  20 

»*yy 


♦  20  D 
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[-  y(. 


cos  20 


■^»5yy 


sin  20  / 


.yyyj] 


<3.17> 


with  / _ and  / _ given  by 

.xyy  »yyy 


^.*yy 


f.yyv 


+i  [*  <x  -  a)2  -  y2 

4nCi  -  v>  [  <c_  .  a>2  +  _*>2 

+2y  [  X  ~  3 

4n<i  -  v)  <c_  .  a>a  +  ^2 


(x  4  a>2  -  y2 
<(x  ■*.  a>2  +  yV 

(3.18) 

x  +  a 

«x  +  a)2  4  y2)2 

<3.19> 


The  functions  /  (x,y,0)  through  /  (x,y,0)  (Chapter 

lx  3y 

II)  can  be  determined  from  the  form  of  the  displacements  and 
stresses: 


u*  *  D. 

“y  “  fa<***>  D. 

o  m  f  <x,y,0>  D 

xx  3x  a 

o'  *  /,  (x,y,0)  D 

yy  4x  a 

T  «  /  (x, y,0)  D 

xy  8x  a 


4  fiy(X,y,0>  Dn 

(3.20) 

♦  f  (x,y,0)  D 

2y  n 

(3.21) 

4  f3y(x,y^>  Dn 

(3.22) 

4  f4y(x,y,0>  Dn 

(3.23) 

4  f  (x,y,0>  D 

By  n 

(3.24) 

The  required  equations  now  exist  for  computer  imple¬ 


mentation  of  the  Displacement  Discontinuity  Method  for  an 
isotropic  plate  in  a  state  of  plane  strain  (program  TVODD  in 
Appendix  C>. 


For  the  case  that  9  =  0°  (along  the  positive  x-axis>,  the 
solution  becomes 


Program  TVODD  was  used  to  solve  this  doubly- symmetric 
problem  using  as  few  as  three,  and  as  many  as  192,  elements 
per  quarter  circle.  Figure  3.3  indicates  that  even  for  the 
case  of  three  elements,  the  solution  obtained  using  the  DDM 
is  extremely  accurate. 

C.  Fracture  Mechanics  Applications 

Since  the  purpose  of  this  study  was  to  validate  use  of 
the  DDM  for  fracture  mechanics  applications,  the  following 
problems  were  solved  using  program  TWODD:  (i)  offset  paral¬ 
lel  cracks,  (ii>  slanted  embedded  crack  (45*  angle),  and 
(iii)  single  edge  crack  near  a  hole. 

1.  Offset  Parallel  Cracks.  The  first  fracture  mechan¬ 
ics  problem  analyzed  was  the  case  of  two  offset  cracks  par¬ 
allel  to  each  other  (Figure  3.4).  The  solution  to  this 
problem  for  an  infinite  plate  was  determined  by  Rooke  and 

Cartwright  (19)  to  be  approximately 

1/2 

Location  “A":  K  ■  1.77  pel  in 
x 

Location  “B“:  -  1.86  pal  in1'2 

with  an  applied  stress  of  1  pel,  and  dimension  "a"  equal  to 
2  inches. 
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Figure  3.4.  Offsai  Paraltol  Cracks  in  an  Infinite  Plata 

(Isotropic  > 


For  the  same  configuration,  the  DDM  yields  the  values 

in  Table  3.1,  indicating  good  convergence. 


Table  3.1.  Offset  Parallel  Cracks  <Isotropic) 


Elements 

K  (A) 

X 

K  (B) 

X 

12 

1.879 

2.023 

48 

1.787 

1 .918 

192 

1.756 

1.886 

2.  Slanted  Embedded  Crack  (45*  Angle).  The  second 
fracture  mechanics  problem  analyzed  was  the  case  of  a  crack 
at  a  45-degree  angle  to  the  direction  of  the  applied  load  in 
a  finite  panel  (Figure  3.5).  Rajiyah  and  Atluri  (16) 
describe  a  boundary  element  alternating  method  which  they 
used  to  obtain  the  solution  to  this  problem: 

K  m  1.378  psi  in1'2 
ii 

when  the  applied  stress  is  1  psi.  For  an  applied  load  of 
1000  psi,  the  solution  would  therefore  be 

K  -  1378  psi  in 
ii 

For  this  1000  psi  load,  the  DDM  yields  (with  192  half¬ 
crack  elements)  the  values  in  Table  3.2,  where  "Length**  is 
the  size  of  each  element  along  the  outer  boundary  of  the 
part  (in  inches).  Again,  good  convergence  is  demonstrated 
—  this  time  for  a  finite  plate. 
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Table  3.2.  Slanted  Embedded  Crack  {Isotropic) 


Length 

K 

z  z 

%  Error 

0.50 

1382 

+0.3 

0.25 

1382 

+0.3 

0.125 

1381 

+0.2 

3.  Single  Edge  Crack  Near  a  Hole.  The  last  fracture 
mechanics  problem  analyzed  was  the  case  of  a  single  crack 
emanating  from  a  hole  at  the  edge  of  a  finite  plate  CFigure 
3.6>.  The  solution  to  this  problem,  according  to  Rajiyah  and 
Atluri  <16),  is 

1/2 

K  -  3.493  psi  in 

i 

under  an  applied  stress  of  1  psi.  Since  this  is  a  Mode  I 
problem,  »  0  psi  in1X2. 

For  a  1  psi  load,  the  DDM  yields  Cwith  96  half-crack 
elements)  the  values  in  Table  3.3,  where  “Length**,  again,  is 
the  size  of  each  element  along  the  outer  boundary  of  the 
plate  <in  inches).  Notice  that  convergence  is  slower  for 
this  more  complicated  boundary.  However,  the  trend  of  mono- 
tonically  approaching  the  solution  as  the  number  of  elements 
increases  is  still  evident. 
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Figure  34).  Single  Edge  Crack  Near  a  Bole  Geotropic) 


Table  3.3.  Single  Edge  Crack  Near  a  Hole  (Isotropic) 


Length 

K 

X 

K  Error 

K 

1 1 

1  .00 

3.785 

+8 .4 

-0.015 

0.50 

3.707 

+6 . 1 

-0.013 

0.^5 

3.663 

♦4.9 

-0.011 

0.125 

3.635 

♦4.1 

-0.009 
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IV.  Orthotropic  Materials 


The  use  of  laminated  composites  in  the  construction  of 
aircraft  components  has  led  to  complications  in  problems  of 
fracture  mechanics.  The  theoretical  solutions  for  isotropic 
materials,  by  themselves,  will  no  longer  suffice  for  deter¬ 
mining  stresses  and  stress  intensity  factors  in  modern 
aircraft  parts.  The  isotropic  DDM,  however,  may  be  modified 
to  include  specially  orthotropic  materials.  That  is 
accomplished  in  this  chapter. 


A.  Theoretical  Development 

Reading  through  Grouch  and  Starfield  <7>,  I  noticed 
that  for  an  isotropic  material  the  form  of  the  displacement 
solution  required  by  the  DDM  could  be  determined  from  the 
solution  to  Kelvin's  problem,  which  involves  a  concentrated 
load  applied  at  a  point  in  an  infinite  solid.  If  the  load 
is  distributed  uniformly  along  the  x-axis  between  x  *  -  a 
and  x  =  +  a  (Figure  4.1),  integration  yields  the  solution  to 
a  constant  traction  applied  to  a  line  crack  <7:48): 


u 

X 


U 

y 


P 

t<3-4 *>/  +  y/vJ  +  -T  t-y/  1 
»y 


p 

l-y/,1  +  2§  t<3-4 y>/  +  y ■/  1 


<4 

<4. 


where  /<x,y>  is  the  same  function 


that  in  Chapter  III  for 


Figure  4.1.  Constant  Traction  Applied  to  a  Line  Grade 


the  isotropic  DDM. 

Differentiating  these  displacements  with  respect  to  y 

only. 


u  "  K4  “  4v>/  ♦  /  1  ♦  3*  I  -/  "  Yf  J 

x,y  20  ,y  ,yy  2Q  ,x  ,*y 


P  P 

u  -  t-/  -  y/  1  IC4-4v>/  4-  /  J 

y.y  2G  J  ,x  TJ,xy  20  J  .y  J  ,yy 


These  may  be  written  in  the  form 


u  -PEA/ 

♦  By/  1  ♦  P  c  C/ 

4*  (Dy/  3 

C4.S> 

x.y  *  ,y 

.yy  y  .x 

.xy 

u  -PIE/ 

y.y  *  .x 

(Fy/  1  ♦  P  [  O/ 
.xy  y  '  ,y 

♦  (Hy/  1 
.yy 

<4.6> 

where  A,  (B,  1H  are  constants  which  depend  only  upon 
material  properties. 


Comparing  to  the  solution  presented  for  a  constant 
displacement  discontinuity  <7:81), 


u  *  D  [2<1  -  v'if  -  yf  ]  +  D  I  -(1  -  2 v>/  -  yf  1 

*  *  ’V  .*x  y  ,x  ,xy 

<4 .7) 

u  ■  D  [<1  -  2v>/  -  yf  3  +  D  [2C1  -  v>/  -  yf  1 

y  *  *X  -xy  y  ,y 

<4.8> 

and  realizing  that  /  =  -/  ,  it  is  seen  that  the  forms 

#XX  iW  ■■ 

are  identical.  Therefore,  in  developing  the  displacement 
discontinuity  solution  for  an  orthotropic  material,  I  used 
the  solution  to  a  constant  traction  along  a  line  crack  as  a 
starting  point  <7:202): 
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<4.9) 


<4.10) 


where 
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\<x  y-r>  ■  jrVfw  -  w] 
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XX 
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T 
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In  r  Cy  > 

2  \ 


C4.11> 


C4.12) 

C4.13> 

C4.14> 

C4.15> 

C4.16> 

C4.17> 

C4.18> 

C4.19> 

C4.20> 

C4.21> 


and  0  is  the  angle  measured  counterclockwise  from  the  global 
x-axis  to  the  line  crack's  local  x-axis.  By  differentiating 
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u  and  u  with  respect-  to  y  only,  we  can  obtain  the  f unc- 

x  y 

tions  required  for  the  DDM  when  the  material  is  orthotropic. 

If  the  form  of  the  solution  for  a  constant  displacement 
discontinuity  along  |  x  |  <  a  is  assumed  to  be 

U*  “  [w  "  @2<jv]  +  ®[w  -  V'V]} 

+  Dy{c[ln  +  D  £ln  r1<^2>/r2<?'2)J|  C4.22) 

“y  -  Dx{E  [w  -  V'.5]  *  -  w]} 

+  Dy|cj^ln  Pi<yi>/Pa<yj> J  +  (H^ln  r1c>'2>/r2<*'2>]}  C4.23) 

then  the  coefficients  A,  B,  ...»  tH  may  be  determined  by 

requiring  u  and  u  satisfy  <i>  the  appropriate  boundary 
x  y 

conditions  and  <ii)  the  equations  of  equilibrium.  Coordi¬ 
nate  transformation  can  then  be  used  to  represent  u  and  u 

x  y 

in  terms  of  the  local  discontinuities  D  and  D  . 

«  y 

The  solution  to  constant  displacement  discontinuities 
D_  and  D_  over  a  line  segment  of  length  2a  Csee  Figure  3.1) 
for  an  ortho  tropic  plate  in  a  state  of  generalized  plane 
stress  is  thus  given  by  Csee  Appendix  A) 
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u 


♦  o  D_ 

X 


+  Q  D. 


-  Q  D_ 
y 


+  Q  D_ 

y 


-  V>v]  -  pi(eI‘>'2>  -  >]]<=“  e 


J* 


In 


-  ln^ri<^>/r2<y2>|  |sin  0 


>]; 


^(W  -  w]  -  -  ®2<’'2>]]sln 

1"(ri<yi>/r2<,'l>)  "  ln(r,<>'2>/r2<>'2:>]jOOS  0 


0 


<4.24> 


u 


y  x 

Y  Y 
1  2 


■  QD. 

X 


12 


r'(ri<?'l>/r2<?'l>]  "  ln(r1<^2>/r2<?'2>)jC:0®  ^ 

-(VJ\>  ‘  W)  “  ~ (W  "  W]  ®ln  * 
1  ^2  J 

r'(rt<ri>/r2<>'l>)  "  ln(*‘1C»'2>/r2<>'2>)|®in  * 


-  Q  D_ 

y 


D_ 

y 


—  fe/V  -  VVl  '  —  (VV  -  02<,'2>  I  lcos  * 

L>V  ^2 


•)]• 


<4.25> 


where 


©  <y  >  ■  arctan 

i  t 


©  <y  >  ■  arctan 

2  t 


r  y/<ytAt>  ] 

[  x  -  a  +  iBJ/A.  J 


<4.26> 


r  ?/c»\v  ] 

iBvy/Av  J 


C4.27> 


*  Cy  >  -  [A 
it  l  t 


(x  -  a)  +  B  (x  -  a)y  ♦  C 


1/2 


<4.28) 
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\>  -  [v 


r  Cy  )  ■  I A  (x  a)2  +  B(x  +  a)y  C  1 1/2 

2  V  I  t  v. 


=/]: 


*■  (j"2  cos2  £  ■*■  sin2  0)  /  y2 
B  -  (1  -  y2'}  sin  20  /  r2 

V  t  v 

C  ■  (.y2  sin2  £  ♦  cos2  0)  /  y2 

L  V  V 


Q  -  — 


12 


2rr  q  y  -  q  y 
1  2  ^2  1 


4  2 

c  c  y  +  [c  (c  +  2c)-cc]y+cc 

11  00  V  12  12  S»  11  22  v  22  60 


(c  y  -  c  )/(c  +  c  > 

11  i  66  12  66 


and  c  are  the  orthotropic  moduli: 


v  E 
12  2 


v  E 
21  1 


11  1  “  t>  V 

12  21 


12  1  -  V  V 

12  21 


1  “  V  V 

12  21 


22  1  L>  V 

12  21 


C  "6 
66  12 


(4.29) 

(4.30) 

(4.31) 

(4.32) 

(4.33) 

(4.34) 

(4.35) 


(4.36) 


The  stresses  for  the  orthotropic  material  are  determined 
from  the  stress-displacement  relationships: 


a 

m 

c 

U 

+  C 

U 

(4.37) 

XX 

11 

x,x 

12 

y.y 

a 

m 

c 

U 

4  C 

u 

(4.38) 

yy 

12 

x.x 

22 

y.y 

T 

m 

C 

U 

♦  C 

u 

(4.39) 

xy 

66 

x.y 

66 

y.x 

For  a  coordinate  system  x,y  at  am  angle  0  with  the  globad 
(material)  coordinate  system,  the  stress-displacement 
relationships  become: 


a  ■  c  cos  6  u  _  -  c  sin  0  u 

xx  11  x,x  11  x.y 

4-  c  sin  0  u  _  ♦  c  cos  flu.  (4.40) 
12  y.x  12  y,y 
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C4.41> 


a  =  c  cos  0  u  _  -  c  sin  8  u 
yy  12  x.x  12  x,y 

+  c  sin  8  u 

22  y,x 

t  ■  c  sin  8  u  _  +  c  cos  8  u 

xy  ®e  x.x  ®®  x,y 

+  C  COS  8  U 

6®  y.x 


C  COS  8  U  _ 
22  ^  y.y 


c  sin  0  u  _ 
®®  y.y 


<4.42> 


Thus,  for  the  case  of  constant,  displacement,  discontinu¬ 
ities,  the  stresses  may  be  written  (Appendix  A) 


D_  [s  c 

X  1  1  1 


cos  8  ~  Sc  sin  8  cos  0 
2  11 


-Sc  sin  0  -  S  c  sin  0  cos  0 

3  12  4  12 

+  (R  -  R  >(c  +  c  /y  y  >sin  0  cos  0 

31  32  11  1212 


-  CR  -  R  >Cc  sin20  -  c  cos20  / y 

41  42  11  12 

-QD_  jSc  sin  0  cos  0  -  S  c  sin 
y  l  1 1  211 

2 

+  S  c  sin  0  cos  0  +  S  c  cos  0 

3  12  4  12 


.V] 


■0 


.4, 

VK  ~ 

3  1 

R  >(c 

3  2  11 

(R  - 

4  1 

R  >(c 

42  11 

a  m  ( 

yy 

0  D*  | 

2 

S  c  cos  | 

- 

S  c 

3  22 

sin20  -  S 

+ 

CR  - 

3  1 

R  >(c 

32  12 

- 

CR  - 

4  1 

R  >(c  1 

42  12 

- 

Q  D_ 
y 

IS  c  sin 

L  1  12 

+ 

S  c 

3  22 

sin  0  cos 

- 

CR  - 

3  1 

R  >Cc  < 

3  2  12 

+ 

CR  - 

41 

R  ><C 

42  12 

12 

♦  c  fy  y  >* 
12  1#  2 

0  -  S  c  sir 
2  12 

c  sin  0  cc 
4  22 

+  C  /}'}'>£ 
22  12 

2 

sin  0  -  c  c 
2  2 

0  cos  0  -  S 
0  S  c  cc 

4  22 

2 

COS  0  -  C  £ 
2  2 

♦  c  >Y  Y  ' 
2  2  '  1*  2 


1*  2 


‘0 


1*  2 


<4.43> 
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T 

xy 


where 


R 

ll 


1  Q  D_  c  S  sin  #  cos  #  *■  S  cos2fi 

x  sail  2 

-  S  sin  0  cos  B  *  S  sin2# 

3  4 

+  <R  -  R  Xsin2#  +  cos2#  /v  y  > 

31  32  *12 

-KR  -  R  XI  -  \/y  Y  >sin  #  cos  #1 

41  42  12  ^ 

-  Q  D_  c#8  ^S^sin2#  ♦  S^sin  #  cos  #  +  cos2# 

-  S  sin  #  cos  #  -  CR  -  R  XI  -  1/v  v  >sin  #  cos  # 

4  3132  1  2  K 

-  <R  -  R  Xcos2#  +  sin2#  /v  Y  >1  <4.45> 

41  42  1  2  I 


s 

1 


s 

2 


s 


3 


s 


4 


©  <r  >  - 

i  i 


e  > 

2  l 


C4.46> 


<4.47> 


<4.48> 


<4.49> 


y  /  Yi 

A  <x  -  a>*  ♦  B  <S  -  a>y  +  Cy2 

v  t.  i 

y  /  r. 

+■  - - - 

A  <Sc  ♦  a>2  t  B  (x  ♦  a>y  +  Cy* 

V  l  1 


<4.50> 
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R  ■  f©  Cy  >  -  ©  Cy  >] 

2i  ^  1  '  i  2  v  J.x 

Cx  -  a>  /  y 

A  Cx  -  a>2  +  B  (x  -  a>y  +  G  y2 

III 

Cx  +  a>  /  y 

A  Cx  +  a>2  +  Bi  Cx  +  a>y  +  C.  y2 

R  ■  In  fr  <*y  >/r  <y  > 

3l  1  i  2  l 

A  Cx  -  a>  +  -B  y 

_  _ v _ 2  i  _ 

A  Cx  -  a>2  +  B  (x  -  a>y  +  C  y2 

v  i  i 

A  (x  +  a)  +  -B  V 

_ v _ 2  v  _ 

A  (x  +  a)2  +  B  Cx  ♦  a>y  ♦Gy2 

V  i  i 

and 


R 

41 


In 


fr  Cj'  >/r  Cy 

[  1  l  2 


4 


A  Cx 


-B  Cx  -  a)  +  G  y 

2  v _ i* _ 

a>2  +  B.  Cx  -  a>y  +  C.y2 

v  t 

-B  <x  +  a)  ♦Gy 

_ 2  V _ l 

A  Cx  +  a>2  ♦  B  <x  ♦  a>y  ♦C.y2 

i  v  v 


Thus,  the  functions  in  Chapter  II  for  an  orthotropic 
material  are 


C4.51> 


C4.52> 


C4.53> 
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/1K<X, y,J0>  -  Q 


In 


-  «*<>%>]  -  pH'W  -  w)]°°*  ' 

Jjsin  Je 


fr  Cy  >/r  <?"  >1  -  in  fr  Cy  >/r  Cy  > 

l  i  »  2  i  J  l  *  *  2  2 


<4.34  > 


/iyCx,y,J0>  -  Q 


-  Q 


In 


(ri<J'1>/ra<,'i>]  "  lcos 


•>); 


— fvv "  V^l 

l  ■  — fw  -  wl 

y  1  1  1  2  1  J 

2 

^1 

j  sin  J0 


<4.55> 


/a x<S,y,J0>  - 


-  Q  I  — 


-  D_  In  [**1<J'1>/Pa<yi>]  “  1n(*,1<J'a>'*,a<J'a>]|c°*  J 

[“fc^  '  W)  "  '  W)]®1"  ^ 


<4.56> 


/2y<S,y,J0>  -  -  Q 


y  ^ 
i  a 


-  w)  -  '  w]]00* 

lnfr  Cy  >/r  Cy  >1  -  Infr  C y  >/r  Cy  >1  I  sin  \ 

^  i ' i  a 'i J  l13  a  '  a  J  I 


0 


J0 


<4.57> 


/  Cx,y,J0>  “  Q  [s  c  cos2  ^  -  S  c  sin  ^0  cos  J0 
7  [_  1  11  2  11 

-Sc  sin2  ^0  -  S  c  sin  J0  cos  J0 
a  12  4  12 

♦  <R  -  R  ><c  +  c  /y  y  >sin  J0  cos  ^0 
»i  aa  ii  12  '  i*  a 


-  <R  -  R  ><c  sin2  ^0  -  c  cos 

41  42  11  ^  12 


2  j@  'rtr,>] 


<4.58> 
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/  <x,y,J0>  ■  -  Q  Sc  sin  J0  cos  ^  -  Sc  sin2  j0 

3y  [  i  u  i  u  r 

♦  Sc  sin  cos  J0  ♦  S  c  cos2 

3  12  4  12 

-  <R  -  R  Xc  cos2  1B  -  c  sin2  JB  /y  y  > 

31  32  11  ^  12  ^  *1*  2 

+  <R  -  R  Xc  +  c  /y  y  >sin  i0  cos  J0|  <4.59> 

41  42  11  12  *1*2 

/  <x,y,J£?>  *»  Q  Is  c  cos2  J/9  -  S  c  sin  3 (3  cos  3  8 

*4X  *  1  12  2  12  K 

-Sc  sin2  J£  -  S  c  sin  38  cos  3B 

3  22  4  22 

+  <R  -  R  Xc  +  c  /y  y  >sin  3&  cos  '’fi 

31  32  12  22  V  2 

-  <R  -  R  ><c  sin2  -  c  cos2  J0  /y  y  >1  <4.60> 

41  42  12  22  1  2  J 

/  <x,y,30>  ■  -  Q  JS  c  sin  cos  J£  -  S  c  sin2  V 

4y  [  1  12  2  12 

♦  Sc  sin  cos  +  S  c  cos2  iB 

3  22  4  22 

-  <R  -  R  >Cc  cos2  30  -  c  sin2  /y  y  > 

31  32  12  22  1  2 

♦  <R  -  R  Xc  ♦  c  /y  y  )sin  30  cos  30\  *  <4.61> 

41  42  12  22  1  2  J 

•*» 

/  <x,y,J£>  »  Q  c  fs  sin  30  cos  V  +  S  cos2  J0 

"X  ww  I  1  2 

-  S^sin  J0  cos  J0  +  S^sin2  V 

♦  <R  -  R  Xsin2  30  +  cos2  J0  /y  y  > 

31  32  K  'l'2 

+<R  -  R  XI  -  1/y  y  >sin  cos  J0|  <4.62> 

41  42  12  K  J 

/5y<x, ,y,30>  m  -  Q  c##  ^S^sin2  ♦  S^sin  '’/S  cos 

+  S#  cos2  V  -  S^sin  V  cos  ^ 

♦  <R  “  R  XI  -  1  fy  y  >sin  30  go s  30 

31  32  12 

-  <R  -  R  Xaos*  V  sin2  3fi  Y  >1  <*.«»> 

41  42  w  F  ftf2  J 

Recalling  also  from  Chapter  II  the  convention  of  com¬ 
puting  element  self-effects  on  the  negative  side  of  the 
discontinuity,  determination  of  the  values  of  coefficients 
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uBgs  and  VBnn  requires  computing  the  following  limits  with 
x  ■  0: 


lylo_  (W  -  W] 

yio.  [ln(''1<>', 


-rr 


0 


<4.64> 

<4.65> 


Thus, 


u 


■  Q  —  <-rt>  -  — ic-rr>  I  fvD_  cos  x0  -  VD_  sin  x0\ 

k  Jl  ’  y  J 

1  ’v,  K  q.l,  ,  >„■> 

*  —  - 1 —  -  — I  (~n>  I  D_  cos  0  -  D_  sin  0 I 

2n  q  y  -  q  y  \Y  Y  J  **  y 

1  2  2  i  w  a  iJ 

>) 


1  q  r  “  q  y  <■ 

- £_i - i-^fvD_  cos  ^  -  VD_  sin 

2  qr  -  qy  <■  X  y 

2  2  1 

2  ^lD_  cos  x0  -  VD_  sin  ’’/sj 


<4.66> 


and 


u  -  -  Q 

y 


C-rt>  -  £vD_  sin  x0  *  ’’D.  cos 


1  YJ%  K  qal  rv  i  v  i  'I 

■  -  —  -  —  -  —  C-n>  I  D_  sin  0  D_  cos  0\ 

2nqy-qyy  y  **  "  y  ^ 

1  2  2  i  t'  i 


1  Va  "  Vid 


£VD_  sin  x0  +  VD_  cos 


2  V*  -  Vi 
i(DS  *ln  ^  +  'Dy  C°* 


C4.67> 


Therefore,  the  local  displacements  become 
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'’u  “  u_  ■  Lu  cos  u  sin  v B 

s  x  x  y 

-  +  -Ceos2  +  sin2  V0>VD 

2  M 

-  +  ~  1 D  <4 .68) 

2  S 

v  “  v  •  i  —  .  v  ^  .  i  •  v  ^ 

u  ■  u_  ■  -  u  sin  6  +  u  cos  B 
n  y  x  y 

t,  .  2  v  .  2  i,  „i  _ 

■  +  -Csin  B  +  cos  0>  D_ 

2  y 

-  +  -  lD  <4 .69) 

2  n 

Thus,  U’B  ■  '■'"B  ■  +  -,  as  stated  in  Chapter  II. 

ss  nn  2 

The  required  equations  now  exist  for  computer  imple¬ 
mentation  of  the  Displacement  Discontinuity  Method  for  an 
orthotropic  plate  in  a  state  of  generalized  plane  stress 
(program  TWODDO  in  Appendix  D>.  This  method  may  also  be 
used  to  determine  the  solution  for  an  isotropic  plate  in  a 
state  of  generalized  plane  stress.  Since  the  ortho  tlropic 
solution  is  singular  for  a  truly  isotropic  material  CE^  m  Ey 

■  E  implies  q  ■  q  ■  1  and  y  ■  y  ),  it  is  not  possible  to 
12  12 

use  E  *  E  *  E  directly.  However,  if  the  values  of  E  and 

x  y  x 

E^  are  within  a  few  hundreths  of  a  percent  of  one  another, 
an  isotropic  material  will  be  sufficiently  represented. 

B.  Stress  Analysis  Application 

To  verify  the  stress  solutions  used  in  program  TWODDO, 
the  case  of  a  circular  hole  in  an  infinite  plate  under  uni¬ 
axial  tension  was  again  considered  (F  ire  4.2).  Lekhnitskli 
<13:163-171)  demonstrates  that  the  tangential  stress  along 
the  boundary  of  the  hole  is  given  by 
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<7t  at  0  =  90  (psi) 


where 


E 

t 


Ai  /Li  -  -  CE  /E  >1/2 
12  1  2 

n  -  <E  /G  -  2v  ~  2fj  V  >1/2 
1  12  12  12 


h 


2 

u  cos  G 
2 


+  <1  n>sin  G 


C4.70> 


<4.71> 

C4.72> 


and 


1 


E 

t 


sin4G 

-  + 

E 

l 


sin2G 


2 

cos  G 


4  _ 

cos  G 


E 

2 


<4.73> 


for  a  uniform  traction  p  at.  infinity. 

For  the  case  G  *  90®,  the  maximum  stress  is  determined 

to  be  ■  a  ■  pCl  +  n>.  Then,  assuming  a  boron/epoxy 

laminate  <E  -  30x10®  psi,  E  -  2.7x10®  psi,  G  -  0.65x10® 

psi,  and  u  «  0.21>,  the  theoretical  maximum  stress  at  the 
12 

edge  of  the  hole  is  a  «  823.9  psi.  The  results  from  pro¬ 
gram  TWODDO  were  extrapolated  to  the  edge  of  the  hole  by 
assuming  an  inverse  fourth-order  least  squares  curve  fit 
<1:264>  is  appropriate.  The  results  are  shown  in  Table  4.1. 


Table  4.1.  Tangential  Stress  at  Edge  of  a  Circular 
Hole  <0rthotropic> 


Elements 

a  at  G»90  * 

t 

X  Error 

3 

304.6 

-63.0 

12 

400.9 

-51 .3 

48 

482.6 

-41 .4 

96 

523.8 

-36.4 

800 

690.0 

-16.3 
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The  extrapolated  stress  at  the  boundary  converges  very 
slowly  due  to  the  presence  of  stress  singularities  at  the 
ends  of  each  element.  The  difference  in  stress  values  a 
short  distance  from  the  boundary,  however,  is  less  than  2% 
when  the  3-element  case  is  compared  to  the  800-element  case. 
Thus,  the  stress  solution  is  accurate  except  very  close  to  a 
boundary  element. 

G.  Fracture  Mechanics  Applications 

Two  of  the  fracture  mechanics  problems  analyzed  for  an 
isotropic  material  were  considered  again  for  an  orthotropic 
material;  (i>  slanted  embedded  crack  (45*  angle),  and  (ii) 
single  edge  crack  near  a  hole.  Analyzed  first  was  the  case 
of  symmetric  edge  cracks  emanating  from  a  hole. 

1.  Symmetric  Edge  Cracks  Emanating  from  a  Hole.  This 
first  fracture  mechanics  problem  involves  cracks  extending 
from  two  sides  of  a  circular  hole  in  an  orthotropic  panel 
(Figure  4.3).  The  solution  to  this  problem  is  presented  by 
Tan  and  Bigelow  (23)  as 

K  =  S  VrTa  F  (4.74) 

z 

where  F  is  a  correction  factor  based  upon  the  geometry  and 
the  material  properties. 

The  case  analyzed  was  a  10*1  graphite/epoxy  laminate 
having  the  following  properties  and  configuration: 
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SYMMETRIC  HOLE  CRACKS  -  ORTHOTROPIC 

3400  q 


3200 

u) 

Cl 


3000 


2800 


2600 


2400 


Border  Element  Size: 
0.125",  0.25",  0.5",  1.0' 


I 


-  11.72  GPa  -  1.70xl0C  psi 
£y  =  144.8  GPa  =  21.0x10*  psi 
G  *9.65  GPa  ■  1.40x10*  psi 

xy 

v  =  0.017  S  -  1000  psi 

xy 

a  =  1.60  in  2a/W  =  0.20  2R/W  =  0.125  H/W  =  2.00 

For  these  conditions.  Tan  and  Bigelow  (23>  predict  a  stress 
intensity  factor  of  ■  2394  psi  in1^2  <F  ■  1.0679>. 

Program  TWODDO  produced  the  results  in  Table  4.2  for  96 
half** crack  elements. 


Table  4.2.  Symmetric  Edge  Cracks  Emanating  from  a 
Hole  <Orthotrop i c > 


Le  ngth 

K 

z 

X  Error 

1  .00 

2327 

-2.8 

0  .56 

2355 

-3.8 

0 .25 

2285 

-4.6 

0 . 125 

2274 

-5.0 

2.  Slanted  Embedded  Crack  (45*  Angle).  The  geometry 
for  this  problem  is  similar  to  that  for  the  isotropic  case 
(Figure  4.4>.  The  material  properties  used,  however,  were 


E 

E 


■  30xl0S  psi 

X 

G  =  0.65x10' 

xy 

•  2.7xlOe  psi 
y 

U  =  0.21 

xy 

which  corresponds  to  a  boron/epoxy  laminate.  Using  a  finite 
element  alternating  method,  Chen  and  Atluri  (4>  determined 
the  Mode  II  stress  intensity  factor  for  the  loading  shown  in 
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SLANTED  EMBEDDED  CRACK  -  ORTHOTROPIC 

<C"'  1.700  q 


a 


Figure  4.4.  Slanted  Embedded  Crack  COrthotropic> 


I 

I 
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Figure  4.4  to  be  ■  1.1097  pel  The  results 

obtained  f rom  program  TWODDO  are  compared  to  this  in  Table 
4.3  for  192  half-crack  elements. 


Table  4.3.  Slanted  Embedded  Crack  (Orthotropic) 


Length 

K 

XX 

X  Error 

0.50 

1  .116 

+0.6 

0.25 

1  .116 

+0.6 

0 . 125 

1  . 116 

+0.6 

3.  Single  Edge  Crack  Near  a  Hole.  Here,  the  same 
geometry  was  used  for  the  orthotropic  case  as  was  used  for 
the  isotropic  case  (Figure  4.5).  The  material  considered  was 
a  C0°1  layup  of  graphite/epoxy  with  the  following  properties: 

E  -  20.5x10*  psi  G  -  0.752x10*  psi 

x  xy 

E  -  1.37x10*  psi  v  -  0.31 

y  xy 

For  this  problem,  Chen  and  Atluri  (4)  obtained  =  3.2349 

i/2 

psi  m  with  a  1  psi  load.  For  an  applied  stress  of  1000 
psi, 

1/4 

K  <■  3235  psi  in  <4.75> 

x 

The  results  from  program  TWODDO  presented  in  Table  4.4  indi¬ 
cate  slower  convergence  than  for  the  more  compact  problems 
analyzed  earlier.  However,  the  trend  of  overestimating  the 


2 


stress  intensity  factor  is  evident  by  the  mo  no  tonic  conver 
gene*  of  the  results. 


Table  4.4.  Single  Edge  Crack  Near  a  Hole  <Orthotropic) 


V.  Conclusions  and  Recommendations 

A.  Conclusions 

As  demonstrated  by  the  result*  in  this  study,  the  Dis¬ 
placement  Discontinuity  Method  may  be  used  as  a  tool  for 
both  stress  analysis  and  fracture  mechanics  applications. 

Its  rapid  convergence  for  problems  involving  infinite 
domains  allows  for  an  accurate  solution  with  only  a  few 
elements,  am  demonstrated  by  both  streets  analysis  applica¬ 
tions  and  by  the  came  of  parallel  offset  cracks  in  an  iso¬ 
tropic  material.  The  method's  accuracy  in  problems  with 
boundaries  composed  of  straight  line  segments  is  adso 
unquestioned.  (See  the  slanted  embedded  crack  results.)  On 
the  other  hand,  for  fracture  mechanics  problems  containing 
curved  boundaries  the  DDM  showed  much  slower  convergence 
than  it  did  in  the  other  cases  discussed  above. 

The  Displacement  Discontinuity  Method  has  an  advantage 
over  finite  element  techniques  due  to  the  fact  that  only  the 
boundary  of  the  problem  must  be  modeled  versus  modeling  the 
entire  domain.  This,  coupled  with  its  rapid  convergence, 
allows  the  DDM  to  obtain  accurate  solutions  using  fewer 
elements  than  finite  element  methods.  Mesh  refinement  is 
also  accomplished  more  easily  with  the  DDM  because  of  the 
reduction  in  order  of  the  problem  from  a  two-dimensional 
domain  to  a  one-dimensional  boundary. 
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The  main  disadvantage  of  the  Displacement  Discontinuity 
Method  is  that  the  matrix  equation  developed  by  the  tech¬ 
nique  contains  a  fully- populated  matrix.  Finite  element 
techniques  are  able  to  use  banded  matrix  operations  to 
increase  their  efficiency;  the  DDM  cannot. 

B.  Recommendations 

This  thesis  dealt  mostly  with  development  and  valida¬ 
tion  of  a  plane  stress  displacement  discontinuity  method. 
Further  research  should  be  done  to  characterize  the  tech¬ 
nique.  This  characterization  could  include  <i>  determining 
the  range  over  which  a  least  squares  fit  is  appropriate  <as 
a  function  of  distance  from  the  crack  tip),  Cii>  obtaining  a 
rough  estimate  of  the  number  of  elements  required  to  obtain 
an  accurate  solution  for  a  given  problem,  and  Ciii)  deter¬ 
mining  the  effect  of  refining  the  mesh  for  a  half-crack  in 
another  manner  Ce.g.  the  method  used  by  Harris  <11  >  for  the 
plane  strain  DDM). 

Another  recommendation  is  to  improve  the  efficiency  of 
Program  TWODDO.  This  could  be  accomplished  by  using  a 
matrix  equation  solution  routine  other  than  Gaussian  elimi¬ 
nation,  which  is  of  order  N*.  Since  the  program  is  modular, 
incorporating  a  different  matrix- solving  algorithm  should 
not  be  difficult. 

Also,  the  current  Displacement  Discontinuity  Method  is 
limited  to  two-dimensional  analyses.  Extension  of  the  DDM 
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to  three  dimensions  may  be  possible  if  a  3- Dimensional 
Fictitious  Stress  Method  could  be  developed. 


1.  Int.roduct.ion 


This  appendix  demonstrates  that  the  orthotropic 
solution  used  in  the  body  of  this  thesis  satisfies  <i>  the 
boundary  conditions  appropriate  4  or  a  displacement 
discontinuity  and  Cii)  the  equations  of  equilibrium  for  a 
panel  subjected  to  a  state  of  generalized  plane  stress. 

The  first  of  these  requirements  (satisfying  the 
boundary  conditions  for  a  displacement  discontinuity >  can  be 
written  in  the  following  form: 


lim  11m 

ii  —  it  ait 

y-»0_  i  y->0+  i  i 


'D. 


)x|<a 


<A.1> 


where  i  »  x,  y.  This  represents  a  discontinuity  in  the 
displacement  field  along  a  line  crack  lying  on  the  x-axis 
between  x  «  -  a  and  x  ■  ♦  a.  This  discontinuity  can  be 
either  a  shear  discontinuity  <i  ■  x>  or  a  normal 
discontinuity  <1  ■  y>. 

The  second  condition  requires  that  the  stresses  satisfy 
the  equations  of  equilibrium: 


a  ♦  t  ■  O  <  x- direction  >  (AJ2) 

kk^c  xy.y 

a  *  r  s  o  <  y- direct  ion  >  (A.3> 

yy.y  *y»* 
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These  equations  are  appropriate  for  a  thin  panel  in  which 
body  forces  and  through- the— thickness  shear  are  neglected. 

The  orthotropic  solution  presented  is  shown  to  be 
correct  as  follows:  First,  the  solution  for  a  shear 
displacement  discontinuity  is  demonstrated  to  be  correct. 
Next,  the  solution  for  a  normal  displacement  discontinuity 
is  validated.  The  two  results  are  then  combined  to  yield  a 
more  general  solution.  This  solution  is  then  generalized 
further  to  determine  the  displacement  field  for  an  arbitrary 
orientation  of  the  crack  with  respect  to  the  material 
coordinate  system.  Finally,  the  stresses  for  this  general 
solution  are  determined. 


2.  Shear  Displacement  Discontinuity 

Let  u^  represent  displacement  in  the  x- direction  and  let 
u  represent  displacement  in  the  y- direction.  Consider  the 
following  displacement  field: 


u 

X 


u 

y 


O  D 

x 


r  r  y/ri\ 

f  y/ym 

1  arctanl - 1 

-  arctanl — - — II 

L  l*  - 

l*  «•  aJJ 

r  r  y/r* 

r  y/ya*nl 

(arctanl - 1 

-  arctan[ — - — 1 1  > 

L  l*  “  •) 

\n  +  J J 

r  r<x  -  a>a  ♦  Cy/r 

— bL.  D  J  in  - — 

x  {  <x  ♦  a>2  ♦  <y/?'1> 


-  In 


Cx  -  a> 


<x  ♦  a> 


♦  <y/rt>a\  ) 


<A.4> 


(A5> 


A-2 


with 


1 

^1 

^2 

2rt 

V*  • 

Vi 

■  <c 

2 

y  - 

C  )/<C 

li  l 

66  12 

■  <c 

2 

Y  “ 

C  >/<C 

11  2 

• 

• 

i* 

4 

r . 

♦  tc 

<c  ♦  2c 

66  V 

12  12  i 

66 


11  22  *  i  22  66 


and 


D  ■  constant 

X 


Than,  on  tha  x-axis 


u  -  Q  D 

X  X 


lint  arctan 
y*»0 


{>. 

-  —  lim  [arctanf  *  — 1  - 

rx  y+o  [  l*  "  */ 


arctan 


arctan 


£*)]} 


U 


<A.6> 

<A.7> 


Tha  limits  da  panel  upon  tha  valua  of  x  and  upon  from  which 
diraction  y  ■  0  is  approached.  The  limits  of  the  arctangent 
function  are 
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This  means  u  ■  0  on  |x|  >  a.  If  |x|  <  a. 
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The  displacements  and  must  also  satisfy  the 
equilibrium  equations: 


a  ♦  t  •  0 

xx,x  xy,y 

O'  ♦  T  =0 

yy#y  *y.x 


<  x- direction  > 
(  y- direction  > 
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which,  when  combined  with  the  stress** strain  relationships. 
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may  be  written 
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Carrying  out  the  differentiation  in  the  x~  equilibrium 
equation. 
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then  each  term  in  Eq  CA.17)  is  zero  and  X- direction  equilib¬ 


rium  im  nti»fi«d.  Using  the  definitions  of  and  it 


is  seen  that 
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If  Eq  (A.22)  holds,  then  x-direction  equilibrium  is  satisfied. 
The  proof  is  given  ir  Appendix  B  CTheorem  1>. 


Carrying  out  the  differentiation  for  the  y- equilibrium 
equation. 
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because  each  berm  in  the  equilibrium  equation  CEq  (A.23»  is 
than  identically  zero. 

Recalling  the  definitions  of  and  q^,  and  realizing 
that  q^  e  i. 
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Thus,  y- direction  equilibrium  is  satisfied  if 
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which  is  the  characteristic  equation. 


3.  Normal  Displacement  Discontinuity 

Again,  let  u^  represent  displacement  in  the  x- direction 
and  let  uy  represent  displacement  in  the  y-  direction.  Now 
consider  a  displacement  field  defined  by 
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with  0  ■  constant,  and  all  other  variables  defined  as  in 

y 

the  shear  displacement  discontinuity  solution. 

Following  the  procedure  above,  for  Ixl  >  a,  u  ■  u  ■  0. 
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The  limits  depend  upon  from  which  direction  y  ■  0  1* 
approached. 
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Thus,  and  satisfy  the  displacement  discontinuity 
boundary  conditions 
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a,  y  -  0 

<A.33> 

u  -  u  ■  0  on  x  > 

y  y  11 

a,  y  ■  0 

<A.34> 

The  displacements  u^  and  must  also  satisfy  the 

equilibrium  equations: 

cu  ♦  c  u  +<c  ♦  c  >  u  « 

11  x,xx  ••  x,yy  12  ••  y,xy 

0  <  x  > 

CA.1S> 

CU  +  C  U  ♦  <c  ♦  c  >  u  ■ 

22  y,yy  ••  y,xx  12  ••  x,xy 

0  <  y  > 

CA.16> 

Carrying  out  the  differentiation  in  the  x~ equilibrium 
equation  (Gq  <A.13», 


c  u  ♦  c  u  ♦  <c 

U  x,xx  ••  *,yy  12 


♦  c  >  u 

••  y,xy 


[<«-•>*-  Cy/r  >2 
- 1 — 

IOC  -  a>*  ♦  Cy/rt>* 3a 

<x  ♦  a>2  -  <y/y  >2  <x  -  a>2  -  <y/r  >* 

♦  - - -  e  - - - 

Kx  ♦  a>*  ♦  C y/y^l*  «x  -  a>*  ♦  <y/r#>*3* 

<x  ♦  a>*  -  <y/y  >*  ' 
s 

Kx  ♦  a>2  ♦  <y/y  >*1* 


♦ 


c  Q  D 

*•  y 


<x  -  a>2  -  <y/yi>*J/>'2 
Kx  -  a>2  «•  <y /y  >*1* 
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t 


»>2  -  <y/rt>*i/r*  c<«  “  •>*  -  <y'r4>a3/y* 


»>2  *  Cy /r>2i2 


<c  ♦  c  >  Q  D  — 
12  ••  y  Yx 


♦  <c  ♦  c  >  Q  D 
12  ••  y 


C<x  -  a>2  ♦  <y/y#>2 J2 


[<x  ♦  a>2  -  <y/yj>2l/y*l 

t<x  ♦  a>2  ♦  <y/y  >V  J 
2  ■■ 


ICx  -  a>  -  <y/rt>  j/y1 
l<x  -  a>2  ♦  (y/y^2]* 


[<x  ♦  •>*  -  <y/r1>2J/rl 

C<x  ♦  a>2  ♦  Cy/y^2!2  J 


tCx  -  a>2  -  <y/y  >*3/y 
2  2 


l<x  -  a>2  ♦  <y/yt>2!2 


♦  a>2  -  Cy/y^^l/yl 

*  ♦  a>*  ♦  <y/y  >212  J 


|C<x  -  a> 


c  +  c  > 

12  •• 


-  !st .  !i<c  *  c  >] 

[Kx  -  a>  <y/y  t>  3  J  t  Yx  Yx  J 

,  <x  *  £,  -  l(c  _  !±±  _  !^0  *  c  >1 

<i  -  »>2  -  <y/y,>* 1  f  <=..  <«.  1 

-  [kx  -  .>*  *  <y/y,>vjt°“  r\  r\  **  ) 

,  '***:*»•*-]{.  . !» - .c  >) 

r  (ICx  ♦  •>*  ♦  Cy/r^vjl  “  r\  y*  “  *•  J 


C  ♦  C  > 
12  •« 


C  +  C  > 

12  •• 


<A.39> 
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Note  that  each  term  is  identically  zero  if 


c  -  c  /y*  -  q  <c  ♦  c  >/y2  ■  0 
11  22  1  1  13  22  *  1 


auid 


c  -  c  /y2  -  q  Cc  ♦  c  >/y*  ■  0 

11  22  2  T  13  8©  2 


Recalling  the  definitions  of  q^  and  q^, 


c  -  c 
n 


•« 


/y2  -  q  <c  ♦  c  >/y2 

f  i  ^1  12  «•  'i 


c  y  -  c 
u'  1  •© 


2 

C  y  -  C  C  +  C 
hi  22  12  22 


c  ♦  c 

12  •* 


-  0 


and 


a  - 
11 


o  /y2  -  qCc  ♦  e  )/y*  ■» 
f3  ^2  i»  ••  r* 


c  y  ~  c 

ir  s  •• 


c  y  -  c  c  ♦  o 
l  la  ••  12  « 


-  0 


c  ♦  c 
12  22 


Thus,  each  term  in  Eq  CA^3)  is  zero,  and  x- direction  equi¬ 
librium  is  satisfied. 

Now,  carrying  out  the  differentiation  for  the  y* 

equilibrium  equation  CEq  <A.16», 


cu  ♦cu  ♦  <c  ♦  c  >  u 

22  y,yy  22  y,xx  12  22  K,xy 


-  c  Q  D  — 
22  y  rt 


2<x  -  aXy *Y  >*Y\ 
[(X  -  a>*  ♦  <y/y  >¥ 


<A.36> 


<a  :arr> 


<A.38> 


CA.39> 
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a <x  -  m >iy*r>*r\  "1 
Kx  -  a>a  ♦  Cy/y^Vj 

q2  '  2Cx  “  aXy/j^J/y* 
rz  ICx  -  a>*  ♦  <y/y2>*l* 

2Cx  -  a>Cy/ya>/ya  ' 

ICx  -  a>a  ♦  Cy/y  >¥ 

2 

qi  2<*  “  *><y/y1> 

ri  t<x  -  a>*  ♦  <y/rx>*f 

2<x  -  a>Cy/yt>  1 

ICx  -  a>a  ♦  Cy/y1>aiaJ 

q,  T  2Cx  -  a><y/y4> 

ra  I  ICx  -  a>*  ♦  C y/r  >¥ 

2Cx  -  a>Cy/ya>  1 
ICx  -  a>a  ♦  Cy/y^VJ 

q  2Cx  -  a>Cy/y  >/y 

y  [  ICx  -  a>a  Cy/yt>ala 

2Cx  ♦  a)(y/j'i>/j'i  2Cx  -  aJCy/^)/^ 

ICx  ♦  a>*  ♦  Cy/ri>ala  ICx  -  a>8  ♦  Cy/ya>*J* 

2Cx  ♦  a>Cy/ya>/ya  1 


2<x  - 


ICx  -  a> 


*><y 'r>'r\  1  f  a 

- Nc  q  - c  ay 

>*  4  Cy/y^Yjl  **  1  M  1  1 


-  Cc  ♦  c  >r 
12  66  1 


2<x  ♦  aXy/y  >/y*  1  f  a 

- Nc  q  -  c  q  y 

[<x  ♦  •>*  ♦  Cy/y  >*1*11  **  1  *•  1  1 

l  ■ 


-  <c  ♦  c  )y 
is  ••  1 


2Cx  - 


CCx  -  a) 


•><y/ya>/y*  1J 

*  ♦  <y/^>Vjl 


"  c6.  V. 


-  Cc  ♦  c 
16  66 


♦  Q  D 


!2Cx  ♦  aXy /Y>'y\  1  [ 

t<x  ♦  *>*  ♦  <y/ya>*l*J\ 


°m\  “  C..  V. 


Cc  ♦  a>rm 
16  66  1 


<A.40> 


Equilibrium  in  ih*  y- direction  i»  aatlafiad  if 


2  2 

c  q  -  c  q  x  -  Cc  +  c  ■  0 

aa^l  66^1r  1  16  66  *  1 


CA.41> 


c  q_  -  c  ay  -  Cc  ♦  c  >y  ■  0 

22  2  66  Te  2  12  M  r2 


CA.42> 


A-18 


since  each  term  in  the  equilibrium  equation  CEq  <A.40>  would 


then  be  zero.  Recalling  the  definitions  of  and  q^. 


2  ,  .  .  I  „ 

c  q  -  c  q  y  -  Cc  ♦  c  ^  >y  ■ 

••  V 1  12  ee  *  l 


Cc  -  c  y  >q 
22  22  i  i 


<c  ♦  c  >y 
12  ••  1 


c  y  -  c 

,  2-  ii*  i  22  ,  .  ..  2 

■  Cc  -  c  *  >  •  Cc  ♦  c  )y 

22  **ri  .  12  22  fX 

C  ♦  C 
12  22 


b-f 

12  22  L 


<c  -  c  ya>Cc  ya  -  c  >  -  Cc  ♦  c  >V* 
22  22  1  11  1  ••  12  22  *  1 


C  C 

11  •• 


y4  ♦  Ic  Cc  +  2c  >  -  c  c  ly*  ♦  c  c 
*1  12  12  ••  11  22  1  22  ee 


c  ♦  C 
12  •• 


CA.43> 


and 


2  ,  .  .*2 

c  q  -  c  q  y  -  Cc  ♦  c  >y  m 

22^2  22^2*2  12  26  €  2 


c  y  -  c 

,  2.  11*  2 

■  Cc  -  c  y  > - 

22  22*2 


Cc  -  c  y  >q_ 

22  22*  2  72 

22  ,  .  v.  a 


Cc  ♦  c  >y 
12  22  *  2 


c  ♦  c 

12  22 


-  Cc  ♦  c  >y 
12  22  *  2 


12 


1 

_  I 

♦  C 

22 


Cc  -  c  ya>Cc  ya  -  c  >  -  Cc  ♦  c  >' 
22  22*2  11*  2  22  12  22 


V*1 


cc  y4  ♦  (c  Cc  ♦  2c  >-cc  ly*  ♦  c  c 
11  22f  2  12  12  22  11  22  *  2  22  22 


C  ♦  C 
12  22 


CA.44> 
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Thus,  y- direction  equilibrium  is  satisfied  if 


c  c  y4<(c  <c  +2c  >  -  c  c  3y24c  c  ■  O  <A.28> 

11  •*'  v  la  la  ••  11  aa  v  aa  •• 


which,  again,  is  the  characteristic  equation. 


4.  Combined  Shear  and  Normal  Displacement  Discontinuities 
Since  the  material  is  assumed  to  be  linearly  elastic, 
the  displacements  resulting  from  one  set  of  boundary 
conditions  may  be  added  to  those  resulting  from  another  set 
to  provide  the  solution  to  a  problem  Involving  both  sets  of 
boundary  conditions.  Therefore,  the  displacements  for 
combined  shear  and  normal  displacement  discontinuities  are 


u  -  ♦  Q  D 
*  x 


{ >  -  — ■»&£=)] 

•  * 


♦  Q  D 


<A.45> 


u  ■ 

y 


f  r<»  -  «>•  -  >‘-r 

l  [<x  +  a>a  ♦  <y/r>*} 

rex  -  a>a  ♦  <y/r9>* 

[<x  a>*  4  <y/y#>*J  J 

\r%  “l  [<*  4  a>*  4 

T<x  -  a>*  4 

[<x  4  a>*  4  Cy/ra>*J  J 


i/’a 


-  in 


A-20 


3.  Coordinate  Transformation  to  Obtain  the  General  Solution 
If  the  crack  ia  not  aligned  with  the  x-axls,  the 
displacements  may  be  determined  by  coordinate 
transformation.  For  a  crack  aligned  with  the  S,  y  system, 
transformation  into  the  x,  y  (material)  system  involves  the 
following  relationships: 


51  cos  0  -  y  sin  0 

CA.47) 

*  sin  0  +  ?  cos  0 

CA.48) 

Uj.  cos  0  -  uy  sin  0 

CA.49) 

u_  sin  0  Uy  oos  0 

<AJSO> 

where  0  is  the  rotation  angle  measured  counterclockwise  from 
the  x-«xls  to  the  5b- axis. 

Carrying  out  the  transformation  on  the  form  of 
the  logarithm  functions  used  in  the  displacement  equations. 


In  lx*  ♦  <y/r>2]1/2  -  In  IA.R2  ♦  B.5??  +  C.yV/2  (A.31) 

A.  -  <r2  com*  0  ♦  sin2  0)  /  y2  <A32> 

V  V  1 

B  -  <1  -  y2>  sin  20  /  y2  (A 

Ci  •  <r2  Sin2  0  ♦  oos2  0)  /  r\  GU54> 


A-21 


leads  'to  using 


In 


A  (Sc  -  a>2  ♦  B  <X  -  a>9  ♦  C  ^l 


A  <X  ♦  a>2  B  <x  ♦  a>9  «■  C  9* 


i/'* 


Although  at  first  it  would  appear  that  the  form  of  the 
arctangent  function  should  be 


r  y/ri  i  r  1  R  *in  ^  +  ?  cos  0  1 

arctanf  - -  I  ■  aratanl  —  -  |  (A.B5) 

x  ^  S  cos  0  -  9  sin  0  *' 

this  would  imply  that  the  crack  lies  along  the  x-axis 
instead  of  along  the  X-axis.  To  account  for  this,  the  form 
of  the  arctangent  function  used  is 


arotan 


l  X  ♦  J 


Thus,  the  following  arctangent  expressions 


appropriate : 


arctan 


\ 


t  <X  -  a> 


aroian 


nv  ] 

*  J 

l  <K  ♦  »>  ♦  ±B.  y/*t  J 


The  limits  as  9  approaches  zero  are  readily  seen  to  be  the 
same  as  before  when  y  approached  zero.  Making  these 
substitutions  yields 


A-22 


.*■<  O 


u  -  ♦  Q  D 

X  X 


and 


l  r*  L  l  <x  -  a>  ♦  iBty/A J 


arcian 


-  arcian 


P/Cy 


f  _^VV_]  ' 

I  <R  ♦  a>  ♦  -B  y/A  J 
v  a  *r  x  * 


(  CS  -  a> 


Y  A  >  *1 

a  2 _  I 

♦  -B  y/A  J 
2  2  2  ^ 


-  arctan 


♦  Q  D 


f  — !^v_  1 1 } 

l  «  ♦  •>  ♦  2B2y/A2  J  J  J 

{fat*  -  a>*  ■*■  B^CR  -  a>y  ♦  CJ2]*'2 

[Ai«  ♦  a>a  ♦  B^CR  ♦  a>y  ♦  C^J 

[A  <R  -  a>2  ♦  B  <R  -  a>y  ♦  C  y2’)1'2  ”1 

- - - - i— [  l  <A36> 

A  <8  ♦  a)  +  B  (S  ♦  a>y  ♦  C  y^J  J 


HA^R  -  a>2  «•  Bt<R  -  a>y  ♦  C^l 
At<R  a>*  ♦  Bt<R  ♦  *>y  ♦  C^J 


IS2 


-  Q  D 


Ia#<r  -  a>  ♦  B#<R  -  a>y  ♦  c^y -I 
A #<R  ♦  a>*  ♦  B  <R  ♦  a>y  ♦  C^J  J 

Jif—f _ !^v_] 

f  *»,v  ]  | 

l  <*  *  •>  *  J  J 

-  —  f  wctanf  _ 

ra  [  l  CR  -  a>  ♦  ;B#y/At  J 

I  —  y/W-  1 1 1 

l  <*  ♦  •>  ♦  J  J  J 


-  arctan 


-  arctan 


<A457> 
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By  making  the  following  definitions. 


©  <r  > 

l  v 


Vn> 


■  arctan 


■  arctan 


f  ] 

l<S  -  a>  ♦  -Bp/A.  J 

f  ) 

(.<*  +  «>♦  ;B.y/*.  J 


r  <r  > 

t  v 


■  (\ 


<S  -  a>  ♦  3 ,<x  -  a>p  ♦  C 


*  <r  >  -  f  A  (S  ♦  a>*  ♦  B 

*  V  ^  V  V 


CS  ♦  a>y  +  C 


f) 

/) 


1X2 


1X2 


<A.38> 


CA  *J9> 


CA.60> 


<A.61> 


the  equations  can  be  written  in  a  more  compact  form: 


u  ■  ♦  Q  D 

X  K 


*  w)] 


♦  Q  D 


u  ■ 

y 


[  r,  (w  *  w)  -  rt  (w 

HS?)-*"!®)] 

[  f;  (w  ‘  w]  -  rt  (w  *  w)j 


-  Q  D 


<A.63> 


Since  the  displacement  discontinuities  of  interest  are 
Dg  and  Dy,  the  above  expressions  are  transformed  using 

■  Dg  cos  0  -  Dy  sin  (J  (A.M) 

D  e  D  sin  0  ♦  D  cos  0  CAJB5> 

y  x  y 
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The  result  is 


u 

X 


-  ♦  Q  D 


I  —  fe  cy  >  -  e  cr  >1 

ra  l  *  1  2  1  J 


-  —  fe  <r  >  -  e  <y  >1 

r,  L  »  a  »  *  J 


-  q  d  —  fe  <r  >  -  e  <r  >1 

y  ra  l  *  1  *  1  J 


'  jr  (w  -  w)  «*"  * 


♦  Q  D. 


H8?HSS]J~ 

•.HWHSS)H 


<xee> 


A-25 


These  equations  are  also  valid  if  the  origin  of  the  R,  y 
system  is  located  at  <c^,c^>  with  respect  to  the  x,  y 
system.  In  that  case,  the  R,  y  coordinates  are  defined  by 

R  ■  (x  -  c  >  cos  0  ♦  <y  -  c  >  sin  0  CA.68> 

x  y 

f  ■  -  (x  -  c  >  sin  0  ♦  <y  -  c  >  cos  0  CA.69> 

x  y 


8.  Stresses  for  the  General  Solution 

The  stresses  are  determined  from  the  st ress** strain 
relationships: 


■ 

o 

c  c  0 

u 

MX 

n  la 

X  ,  X 

o 

.  m 

c  c  O 

* 

U 

•  yy 

la  aa 

y  *y 

T 

0  0  c 

u  ♦  u 

xy  J 

•  • 

■ 

x , y  y,x 

<A.14> 


and  the  coordinate  transformation  relationships: 


_9 

9k 

_9 

*Y 


(cos  0>  -J~  -  (sin  0>  ^ 
<sin  0>  ♦  (cos  0>  ^ 


(A.70> 


<A.71> 


These  are  confined  to  yield 


o  *c  cos  0  u  _  -  c  sin  0  u 

XX  11  X,R  11  ! 


k,9 


4  c  sin  0  u  _  ♦  c  cos  0  u  „  <A.72> 

ia  y,S  la  ^  y,9 
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ct  m  c  cos  8  u  _  -  c  sin  8  u  _ 
yy  12  x,x  12  x,y 


♦  c  sin  fl  u  ,  +  c  cos  £  u  _  <A.73> 

22  y»s  a*  y.9 


t  sc  sin  S  u  c  cos  u  „ 
xy  M  x,x  ee  x,9 


♦  c  cos  8  u  _  -  c  sin  (9  u  _ 
66  ^  y»x  66  y,y 


<A.74> 


With  th«  definitions  of  u  and  u  above,  we  determine  the 

x  y 

following  partial  derivatives: 


m 


^  e<r  > 

dR  2  v 


o  <r> 


d?  a 


e  <y  > 

a  i 


as  ln(r,<J'iVr1<n>) 

ay  ln(r,<^i>/r2<n>) 


■  as  *"(*,,<n>)  -  as  ln(v>\>) 

■  ay  *  ay  “(w) 


These  are  readily  obtained 


e  <y  >  - 

dR  i  fi 


<R  -  • 


?/r  a. 

_  V  V _ 

♦  -b  y/A  >a  ♦  <y/r  a  >* 
a  x  v  v  v 


<A.7S> 


3*  " 


♦  Vy/A.j*  ♦  <y/y.A>a 

a  v  v  v  v 


CA.76> 


•W  - 


<R  -  a>/y,  A. 

I  % 


oz  -  a  ♦  ~b  y/A  >a  ♦  <y/r.A  >a 

i  i'  V  '  x  x 


<A.77> 
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CA.78> 


9 

*9 


e  <r  >  - 

*  '  V 


♦ 


CB  4*  a>/r  A. 

V  V 

Oc  +  •  ♦  -B  y/ A  >*  +  Cf/r  A  >* 
2  i7  v  7  'v  v 


«f  *n(r.<n>) 

<4  ln(rj<n>] 

4  ln(,‘.<n:>) 

4  *n(r/n>] 


A  CR  -  a>  ±B.  p 
_ i _ a  v  _ 

A  CS  -  a>*  ♦  B.<S  -  a>9  ♦  C.9* 

V  V  V 

A.  <8  +  a>  +  iB  9 
_ v _ a  v7 _ 

AO?*  •>*  ♦  B  a  +  «>9  ♦  C  9* 

V  V  V 

-B  <R  -  e>  C  9 

_ 2  V _ v_ _ 

A  CE  -  a>*  B  CE  -  a>9  ♦  C.9* 

V  |  l 

-B.  CS  ♦  a>  ♦  C  9 
_ a  v _ v7 _ 

AO?  *  a>*  ♦  B.<8  *  a>9  ♦  C.9* 

V  V  V 


But.,  from  Thoorom  2  in  Appendix  B, 


CS  ± 


-B  9/A  >*  ♦  C9/y.  A>* 


2  v 


V  V 


A  CR  ±  a>*  ♦  B  CE  ±  a>9  ♦  0.9* 


So  the  derivatives  of  9  and  ©#  can  be  written 


9 

Wt 


e  <r> 

I  V 


9'r. 

A  CS  -  •>*  +  B.  OK  -  *>9  ♦  C  9* 

V  V  V 


9 

9R 


©  Cy  > 
a  v 


Wri 

A  CE  ♦  a>*  ♦  B.  CB  *>9  ♦  C.9* 

V  l  X 


9 

3? 


©  Cy  > 

1  V 


CE  -  a>/y. 

A  CB  -  •)*  ♦  B  CR  -  a>9  ♦  C.9* 

V  X  V 
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CA.79> 


CA  JBQ> 


CA.8i> 


CA.82> 


CB.14> 


<A.®3> 


CA.84> 


CA.8S> 


CA.86> 


_  at  ♦  miy/r 

-=  ©  <y. )  ■  +  - - - — 

^  2  v  A  <Sc  ♦  a>2  ♦  B.  <5c  ♦  a>y  ♦  C.y 

til 

Thus,  we  can  define  the  following  tar  ms: 

r  .  -  -=£[©  <r  >  -  e  <r  >l 

li  aR|_  1  'i  2  I  J 

9/rv 

«  —  ■  -  ■  —  ■  -  ■  — — 

a  cs  -  •>*  ♦  b  <s  -  a>y  «•  c.y2 

V  V  i 

9'r. 

- 

A  <51  ♦  a>2  ♦  B.  «  ♦  a>y  ♦  C.y2 

V  V  I 

r  .  ■  ^fe  c r>  -  ©  <r  >1 

2i  1  S  2  'v  J 

or  -  a >/r 

m  - 

A  < Jt  -  a>2  «■  B.  CR  -  a>y  ♦  C.y2 

V  V  V 

<R  ♦  a)/r. 

^  _ V 

A  O?  ♦  a>2  ♦  B.  Ol  ♦  a>y  ♦  C.f* 

t  V  V 

R  -  In fr  <r>/r  <r  >1 
*i  AS  ^  1  ‘i  2  i  J 

A.  <Jc  -  a>  ♦  ^B.y 
i  2  i* 

■i  ■  ■  — 

A  <R  -  a>2  ♦  B.  <R  -  a>y  C.y2 

i  i  V 

A.  Ot  ♦  a>  ♦  “B.  y 
_ v _ a  v  _ 

A.<*  ♦  a>*  ♦  BC<JI  ♦  *>y  ♦  c^y* 


A  “29 


CA.87> 


CA.88) 


CA.89> 


t  ■  lnfr  Cy.)/r  >1 
4v  ^  i  v  a  v  J 


-b  <s  -  «>  ♦  c.y 

a  v  v 


A  <R  - 

V 


»>*  ♦  b.<r  -  a>y  c.y* 

V  V 


-B.  <R  ♦ 
a  v 


k>  ♦ 


c.y 


A  CR  ♦  a>  ♦  B.  <R 

V  v 


a >9  ♦  G.y2 


<A.90> 


Wow,  with  Eqs  <A.66>  and  <A.67>  we  have 
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♦  Q  D, 


[R..-  0  * Q  d*K'  r“] 


win  ^ 


ein  0  <A.91> 


u  •  Q  D. 

*,Y 


-f-R  *-R  1 

5  r •  11  “J 

[R*““  R~3 


coe  0  -  Q  D_ 


R**‘  R-J 


ein  0 


Q  D 


ooe  0  ♦  Q  D 


(V  R*J 


ein  0  CA.92> 


y.S 
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-  Q  D, 


ft 

<1- 

ft 

q.  1 

—  R  - 

—  R  coe  0  - 

0  D= 

—  R  - 

—  R 

fa  “ 

r,  -J 

R 

r%  » 

r*  iaJ 

ein  0 
CA.93> 


y.5?  r 


£  ■>«[*..-  *-]— *  -  D?(v  R«hn  * 

9^  R«-  J  R~]°~  ^  ^  R«] 


-  Q  Dr 


ein  0 
<A.94> 
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Defining 


%  qi 

S  -  —  R - -  R 

1  r  «  r  « 


—  R 

r*  21 


—  R 


<1  q 

S  ■  —  R - -  R 

s  yx  n  r2  I* 


1R 
''i  21 


Yx  « 


the  partial  derivatives  become 


Q  D_S  cos  0  -  Q  D  S  sin  0 

x  1  y  2 


Q  DCR  -  R  >cos  0 
9  91  92 


Q  D  CR  -  R  >sin  0 

S  91  92 


CA.93> 


u„  y  ■  Q  drs,cos  0  -  Q  DySasin  0  ♦  Q  DyCR^-  R<a>cos  0 


♦  Q  D  CR  -  R  >sin  0  CA.96> 

S  41  42 


u  _  -  -5-  D^<R_  -  R  >cos  0 - —  D  (R  -  R  >sin  0 

%/  S  v  v  8  »4  ee  ^  v  v  R  »4  A 4ft  “ 


y*5  8  •»  *2 

i  a 


y,r2  ?  *i  ** 


-  Q  D?S#co s  0  -  Q  DRS#sin  0 


ca.97> 


D  CR  -  R  >cos  0 - —  D  CR  -  R  >sin  0 

J'i?'a  X  41  42  ^  ytra  9  41  42  ^ 


-  Q  D_S^cos  0  -  Q  DgS4sin  0  CA36) 


Using  the  stress-strain  relationships  CEq  CA.14»  results  in 


A-01 


c  cos  B  u  -  c  sin  B  u  _  +  c  sin  B  u  _  ♦  c  cos  B  u 

44  “  >.  S  44  ^  W  R  ««  V*  3  « 4  « 


ii 


M,X 


11 


ia 


y*a 


ia 
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-  Q  D 
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c  sin  0  cos  fi  -  S  c  sin  ^ 
ii  a  ii 


♦  S  c  sin  0  cos  fl  *  S  e  cos  0 

a  la  4  u 

-  CR  -  R  ><c  cos2  0  -  c  sin2  0  /y  y  > 
ai  sa  ii  ia  r  't'» 

♦  CR  -  R  >Cc  ♦  c  fy  y  >sin  0  cos  pi 

4i  4a  ii  ia  *  i*  a  J 
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4  4k  “  u  C  « 4%  ”  *4  R  SO  r  y  Q  r 


12 
-  Q  D 


E* 


X* 
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X.9 


y.9 


is 


0  -  S  c  sin  0  cos  0 
a  ia 


-Sc  sin  0  -  S  c  sin  0  cos  0 
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si  sa  is  ->■  as  *  i*  a 


-  <"41-  R4a><°i.*ln  *  - 


0  /r 


-  Q  d 


c  sin  0  cos  0  -  S  c  sin  0 
is  a  la 
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a  aa  K  4  aa 


-  CR  -  R  >Cc 
si  sa  is 


*  0  -  c  sin2  0  'y.Y> 

22  2  3 
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4i  42  la  aa  #  r  a  J 
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*y 
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♦  CR  -  R  >Csin  0  ♦ 
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-  Q  Dc 
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CA.101> 


As  can  be  seen  from  the  definitions  of  the  R.  's  in  Eqm 

VJ 

(A.87>  through  <A.90>,  the  stresses  tend  to  zero  as  the 
distance  from  the  origin  tends  to  infinity.  Thus,  the  final 
boundary  condition  of  zero  stress  at  infinity  is  satisfied. 


Appendix  B:  Proofs 


Theorem  1.  Given  the  material  constants  c  ,  c  ,  c  ,  c  , 
-  11  la  aa  ee 

Y  »  Y  and 
*  i  'i 


c  Y  “  c 
li  l  ee 


C  Y  -  C 

112  ee 


c  ♦  c 
12  ee 


c  ♦  c 
12  ee 


then  ■  1. 


Proof: 


Using  the  definitions  of  q^  and  leads  to 


<V»a 


,  2  2  . 
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12  ee 


Define 


62  -  Ic  c  -  c  <c  ♦  c  >1*  -  4  c  c  c2  . 
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Then,  from  the  characteristic  equation: 


4  2 
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n  ee  ia  12  ee  ii  n  22  ee 


we  have 
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Theorem  2.  Given  the  material  const-ante  r  ,  Y  >  the 
-  t  a 

measured  angle  0,  and  the  half- crack  length  a,  and  defining 
the  parameters 

■  <y*  cos2  0  ♦  sin2  0>  / 

B  -  <1  -  y2>  sin  20  /  r* 

V  V  V 

C.  ■  <y2  sin2  ^  ♦  cos2  0">  /  y2 

then 
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2  i  v  v  t 
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t  V  '  V  l 

where  S,  y  are  the  local  coordinates. 

Proof: 

Expanding  the  binomial, 

(52  ±  a  +  -B  y/A.J2  ♦  cy/ty.A.JJ2 
2  i/  i  9  v  v 
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ii  4  v  vv 

Using  the  definition  of  B^: 

B  ■  <1  -  r2>  ein  20  /  y\  »  2  <1  -  sin  0  cos  0  /  r\  <B.10> 
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while  using  the  definition  of  G.  load*  to 
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Appendix  C:  Program  TVODD 


NOTE:  The  information  in  this  appendix  is  taken  directly  from 
Crouch  and  Starfield  (Reference  7)  with  only  slight  modification. 

1.  Input  (Data)  Deck 

Program  TWODD  requires  the  following  six  sets  of  cards: 

Set  1:  Free  FORMAT;  one  card  must  be  given  (may  be  blank). 

Columns  1-80  of  this  card  contain  any  desired  information  to 
identify  the  problem  being  solved. 

Set  2:  FORMAT  (314);  one  card  must  be  given  to  specify  the  following 
control  parameters. 

numbs  -  number  of  straight  line  boundary  segments  (each 
containing  at  least  one  boundary  element)  used  to 
define  boundary  contours. 

Ninos  *  number  of  other  line  segments  (not  on  a  boundary) 
along  which  displacements,  stresses,  and  stress 
intensity  factors  are  to  be  computed. 

1  no  symmetry  conditions  imposed. 
ksym  *  2  x=XSYM  (Card  3)  is  a  line  of  symmetry. 

3  y*YSYM  (Card  3)  is  a  line  of  symmetry. 

4  x*XSYM  and  y*YSYM  (Card  3)  are  lines  of  symmetry. 

Set  3:  FORMAT  (F6.2,E11.4,2F12.4);  one  card  must  be  given  to  define  the 
elastic  constants  and  specify  the  locations  of  lines  of  symmetry 
(if  any). 

PR  *  Poisson's  ratio  (v»). 

E  *  Young's  modulus  (E). 

XSYM  *  location  of  line  of  symmetry  parallel  to  y-axis  (XSYM 
is  ignored  if  KSYM  »  1  or  3  on  card  2). 
ysym  *  location  of  line  of  symmetry  parallel  to  x-axis  (YSYM 
is  ignored  if  KSYM  *  2  or  4  on  card  2). 

Set  4:  FORMAT  (3E11.4);  one  card  must  be  given  to  define  the  initial 
stresses  (if  any)  in  the  region  of  interest. 

PXX  *  <7  at  infinity. 

KM 

PYY  »  a  at  Infinity, 
yy 

PXY  *  r  at  infinity, 
xy 


0-1 


set  5: 


FORMAT  (I4,4F12.4,I4,2E11.4);  NIMBS  cards  must  be  given  to 
define  the  locations  and  boundary  conditions  of  the  boundary 
elements. 


NUM  =  number  of  equally  spaced  boundary  elements  along  a 
straight  line  segment,  all  elements  having  the  same 
boundary  conditions. 

XBEG  =  x -coordinate  of  beginning  of  line  segment. 

YBEG  a  y -coordinate  of  beginning  of  line  segment. 

XEND  =  x -coordinate  of  end  of  line  segment. 

YEND  a  y -coordinate  of  end  of  line  segment. 


KODE 

BVS 

BVN 


'  1  a  and  prescribed. 

2  u*  and  un  prescribed. 

3  us  and  a?  prescribed. 

4  c/%  and  un  or  escribed . 

•  n 

resultant  shear  stress  (<?  )  or  displacement  (u  ). 
resultant  normal  stress  (o^)  or  displacement  (u^). 


Set  6: 


format  (4F12. 4,214);  numos  cards  aust  be  given  to  define 
locations  of  points  inside  the  region  of  interest  vhere 
displacements,  stresses,  and  stress  intensity  factors  are  to  be 
computed. 

XBEG  a  x -coordinate  of  first  point  on  line. 

YBEG  *  y-coordipate  of  first  point  on  line. 

xend  »  x -coordinate  of  last  point  on  line. 

YEND  «  y-coordinate  of  last  point  on  line. 
nu*c>b  *  number  of  equally  spaced  points  between  the  specified 
first  and  last  points. 


KTYPE 


0  compute  displacements  and  stresses  only. 

1  compute  stress  intensity  factors  as  veil. 


Note:  If  KTYPE>1,  (XBEG, YBEG)  must  be  the  location  of 
a  crack  tip,  and  the  line  defined  by  (XBID,YEMD)  aust 
extend  in  the  direction  of  the  crack  in  order  for  the 
computed  stress  intensity  factors  to  have  any  meaning. 


nono 


2.  Program  TWOOP  Listing 

PROGRAM  TWODD 
C 

CCfMCN/Sl/PI,PR,PRl,PR2,C0N,CCNS 

OCm»/S2/SXXS ,  SXXN ,  SYYS ,  SYYN ,  SXYS ,  SXYN,  UXS  ,  UXN ,  UYS  ,  UYN 
CO*t10N/S3/C(  1600, 1600 ) ,B( 1600 ) , D( 1600 ) 

C 

DIMENSION  XM( 800 ) , YM ( 800 ) , A( 800 ) , OOSBET ( 800 ) , SINBET( 800 ) ,KOD( 800 ) 
DIMENSION  TITLE( 20) 

C 

OPEN  (UNIT-5, FILE*' TWODD. DAT', STATUS-  'OLD') 

OPEN  (UNIT-6,FILE*'TNC>OO.LIS ' , STATUS- 'NEW ' ) 

READ  (5,1)  (TIT1£(I),I=1,20) 

WRITE  (6,2)  (TITLE(I), 1-1,20) 

READ  (5,3)  NUMBS, NUMQS, KSYM 
READ  (5,4)  PR,E,XSYM, YSYM 
READ  (5,5)  PXX,PYY,PXY 
WRITE  (6,6)  NUMBS, NUMQS 
GO  TO  (80,85,90,95), KSYM 
80  WRITE  (6,7) 

GO  TO  100 

85  WRITE  (6,8)  XSYM 
GO  TO  100 

90  WRITE  (6,9)  YSYM 
GO  TO  100 

95  WRITE  (6,10)  XSYM, YSYM 

C 

100  CONTINUE 

WRITE  (6,11)  PR,E 
WRITE  (6,12)  PXX,PYY,PXY 
C 

PI-4.*ATAN(1.) 

OON«1./(4.*PI*(1.-PR) ) 

CONS-E/U.+PR) 

PR1-1.-2. *PR 
PR2»2.*(1.-PR) 

DEFINE  LOCATIONS,  SIZES,  ORIENTATIONS  AND  BOUNDARY  CONDITIONS  OF 
BOUNDARY  ELEMENTS. 

NUMBE-0 

DO  110  N-l, NUMBS 

READ  (5,14)  NUM,XBBG,YBBG,XEND,YEND,KODE,BVS,BVN 
XD-(XEND-XBBG)/NUM 
YD-  ( YEND-YBBG )  /NUM 
SW-SQRT ( XD*XD+YD*YD ) 

DO  110  NE-1,NUM 
NUMBE-NUMBE+1 
M-NUMBE 


03 


XM( M ) -XBEG+O . 5* ( 2 . *NE-1 . ) *XD 
YM(M)=YHEG+0 . 5*  ( 2 .  *NE-1 .  )*YD 
A(M)=0.5*SV 
SINBET(M)«YD/SW 
OOSBET(M)-XD/SW 
KOD(M)*KODE 
m«2*M 
MS=MN-1 
B(MS)*BVS 
110  B(MO*BVN 
C 

WRITE  (6,13) 

DO  115  M»1,NUHBE 
SIZE-2. *A(M) 

ANGLE=180. *ATAN2(SINBET(M) ,COSBET(M)  )/PI 

WRITE  (6,15)  M, KOD(M),XM(M),YM(M), SIZE, ANGLE, B(2*M-1),B(2*M) 
115  CONTINUE 
C 

C  ADJUST  STRESS  BOUNDARY  VALUES  TO  AOOOUNT  FOR  INITIAL  STRESSES. 
C 

DO  150  N-1,NUMBE 
NN=2*N 
NS-NN-1 
OOSB-OOSBET(N) 

SINB=SINBET(N) 

SIGS- (PYY-PXX )  *SINB*OOSB+PXY*  (COSB*OOSB-SINB*SINB> 
SIGN-PXX*SINB*SINB-2. *PXY*SINB*C0SB+PYY*O0SB*O06B 
GO  TO  (120,150,130,140),KOD(N) 

120  B(NS)»B(NS)-SIG6 
B(NN)»B(NN)-SIGN 
GO  TO  150 

130  B(NN)»B(NN)-SIGN 
GO  TO  150 

140  B(NS)»B(NS)-SIGS 
150  CONTINUE 
C 

C  COttHJIE  INFLUENCE  COEFFICIENTS  AND  SET  UP  SYSTEM  OF  ALGEBRAIC 
C  EQUATIONS . 

C 

DO  300  I sl, NUMBS 
IN*2*I 
IS-IN-1 
XI*XM( I ) 

YI-YM(I) 

COSBI-COSBET(I) 

SINB1-SINBET(I ) 

KOOE-KOD(I) 

C 

DO  300  J-1,NUMBE 

JN«2*J 

JS-JN-1 
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CALL  INITL 
XJ=XM(J) 

YJ*YM(J) 

COSBJ-COSBET(J) 

SINBJ-SINBET(J) 

AJ=A{ J ) 

CALL  CQEFF(XI,YI,XJ,YJ,AJ,C0SBJ,SINBJ,+1) 

00  TO  (240,210,220,230), KSYM 

210  XJ*2 . *XSYM-XM( J ) 

CALL  CCEFF(XI,YI,XJ,YJ,AJ,COSBJ,-SINBJ,-l) 

GO  TO  240 

220  YJ«2.*YSYM-YM(J) 

CALL  OQBFF(XI,YI,XJ,YJ,AJ,~COSBJ,SINBJ,-1) 

GO  TO  240 

230  XJ=2.*XSYM-XM(J) 

CALL  OQEFF(XI,YI,XJ,YJ,AJ,COSBJ,-SINBJ,-l) 

XJ«XM(J) 

YJ*2 . *YSYM-YM(  J ) 

CALL  GQEFF ( XI ,  YI ,  XJ ,  YJ ,  AJ  , -COSBJ , SINBJ , -1 ) 

XJ*2.*XSYM-XM(J) 

CALL  C0EFF(XI,YI,XJ,YJ,AJ,-G06BJ,-SINBJ,+1) 

240  CONTINUE 

GO  TO  ( 250, 260, 270, 280 ) ,KODB 

250  C(IS,  JS)-(SYYS-SXXS)*SINBI*CQSBI+SXYS*(O0SBI*O06BI-BINBI*SINBI ) 
C(IS,JN)»(SYYN-SXXN)*SINBI*OOSBI+SXYN*(OOSBI*GOSBI-SINBI*SINBI) 
C( IN, JS )«SXXS*SINBI*SINBI-2 . *SXYS*SINBI *OOSBI +SYYS*C08BI *OQSBI 
C(  IN,  JN)*SXXN*SINBI*SINBI-2  .*8XYN*SINBI  *COSBI+SYYN*008BI*G06BI 
GO  TO  300 

260  C(IS,JS)«UX8*OQBBI+UYS*SINBI 
C(  IS, JN) «UXN*C08BI  HJYN*SINBI 
C(IN,JS)«- UXS*SINBI+UYS*COSBI 
C(  IN, JN) — UXN*SINBI +UYN*C08BI 
GO  TO  300 

270  C(IS,J8)«UXS*C06BI+UYS*SINBI 
C(  18,  JN )  -UXN*C08BI+UYN*SINBI 

C(  IN,  JS  )«SXXS*SINBI*SINBI-2 .  *SXYS*SINBI*C0SB1+SYYS*C0SBI*CQSBI 
C(  IN,  JN )  -SXXN*SINBI  *SINBI  -2 .  *SXYN*SINBI  *COSBI  +SYYN*COSBI  *G06BI 
GO  TO  300 

280  C( I 8 , JS ) * ( SYYS-SXXS ) *SINBI *OOSBI +SXYS* ( COSBI *006BI -SINBI *8INBI ) 
C(  IS ,  JN ) « ( SYYN-SXXN )  *SINBI  *008BI +SXYN*  (OOSBI  *OOSBI -SINBI  *SINBI ) 
C( IN, J8 ) «-UXS*SINBI 40X8*00861 
C(IN,JN)  — UXN*SINBI4UYN*C0SBI 


OS 


flOO  o  o  o 


300  CONTINUE 


SOLVE  SYSTEM  OF  ALGEBRAIC  EQUATIONS. 


N=2*NUMBE 
CALL  SOLVE(N) 

COMPUTE  BOUNDARY  DISPLACEMENTS  AND  STRESSES 

WRITE  (6,16) 

DO  600  Ial,NUKBE 
IN»2*I 
IS-IN-1 
XIaXM(I) 

YI«YM(I) 

COSBI  =COSBET ( I ) 

SINBI*SINBET( I ) 

UXNEG-0. 

UYNBG-O. 

SIGXX-PXX 

SIGYY=PYY 

SIGXY-PXY 

DO  570  Jal,NUMBE 

JN»2*J 

JS«JN-1 

CALL  INITL 

XJaXM(J) 

YJaYM(J) 

AJaA( J ) 

COSBJ-OOSBET(J) 

SINBJ «SINBET(J) 

CALL  OOEFF ( XI , YI , XJ , YJ ,  AJ , GOSBJ , SINBJ , +1 ) 

GO  TO  ( 540, 510, 520, 530 ) ,KSYM 

510  XJ=2.*XSYM-XM(J) 

CALL  C0EFF(XI ,YI ,XJ,YJ,AJ, C06BJ, -SINBJ, -1) 
GO  TO  540 

520  YJ»2.*YSYM-YM(J) 

CALL  COEFF(XI,YI,XJ,YJ,AJ, -COS BJ, SINBJ, -1) 
GO  TO  540 

530  XJ*2. *XSYM-XM( J ) 

CALL  C0EFF(XI, YI,XJ,YJ,AJ,COSBJ, -SINBJ, -1) 
XJaXM(J) 

YJ-2.*YSYM-YM(J) 

CALL  C0SOT(  XI,  YI,XJ,YJ,AJ,-C08BJ,  SINBJ, -1) 
XJ«2 .  *XSYM-XM(  J ) 

CALL  CCeFF(XI,YI,XJ,YJ,AJ,-G08BJ, -SINBJ,  +1) 


c 

540  CONTINUE 
C 

UXNBO»UXNBG+UXS*D  ( JS )  HON*D(  JN ) 

UYNEG“UYNEG+UYS*D(JS)+UYN*D(JN) 

SIGXX*SIGXX+SXXS*D( JS)+SXXN*D(  JN) 

SIGYY-SIGYY+SYYS*D( JS )+SYYN*D( JN ) 
SIGXY«SIGXY+SXYS*D(JS)+8XYN*D(JN) 

C 

570  CONTINUE 
C 

USNE0-UXNEG*CO6BI  +UYNEG»8INBI 
UNNBO-UXNEQ*SINBI  +UYNBG*C06BI 
USPOSHJSNEG-D(IS) 

UNPOS-UNNBG-D(IN) 

UXPO6«USPOS*C06BI  HJNP06*SINBI 
UYPOS-USPOS*SINBI-HJNPOS*COBBI 

SIGS*  ( SIGYY-SIGDOC )  *SINBI  •C06BI+SIGXY*  (COSBI  *006BI-SINBI  *SINBI ) 
SIGNsSIGXX*SINBI  *SINBI  -2 .  *8IGXY*SINBI  *006BI  +SIOYY*OOSBI  *C06BI 
C 

WRITE  (6,17)  I,D( IS)  ,USNBG,USP06,D(IN)  ,UNNBG,UE>06,UXNBG,UYNBG, 
1  UXPOS,UYPOB ,310b, SIGN 

C 

600  CONTINUE 
C 

C  COMPUTE  DISPLACEMENTS  AND  STRESSES  AT  SPECIFIED  POINTS  IN  BODY. 

C 

IF  (NUMOS.LE.O)  GO  TO  900 

NPOINT-O 

DO  900  N«1,NUM06 

READ  (5,19)  XBBG,YBBO,XEND,YB©,NUM>B,KTYPE 
IF  (KTYPE.BQ.O)  WRITE  (6,22) 

IF  (KTYPE.EQ.l)  WRITE  (6,18) 

DELX«(XH4D-XBBQ)/NUN? 

DELY- ( YB4D-YBB0) /NUM> 

IF  (NUVB.OT.O)  NH*»«NUM»+1 
IF  (DELX**2-HDBiY**2.EQ.O.)  NUMP-1 
C 

ANGLE»ATAN2  ( YEND-YBEG,  XEMD-XBBO) 

C 

DO  890  NI-1,NUK> 

XP«XBBG+(NI-1)*DKLX 

YP*YBBG+(NI-1)*MLY 

C 

UX-0. 

UY-0. 

SIGXX-PXX 

SIGYY-PYY 

SIGXY-PXY 

C 
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DO  880  J«1,NUMBE 

JN»2*J 

JS=JN-1 

CALL  INITL 

XJ-XM(J) 

YJ»YM(J) 

AJ«A(J) 

C 

IF  (SQRT((XP-XJ)**2+(YP-YJ)**2).LT.2.*AJ)  GO  TO  890 
IF  (KTYPE.E3Q.1)  R*SQRT(  ( XP-XBSG ) * ( XP-XBBG )  +  ( YP-YEEG ) * ( YP-YB8G ) ) 
C 

OOSBJ-QOSBET(J) 

SINBJ*SINBET( J) 

CALL  CQEFF(XP,YP/XJ,YJ,AJ,COSBJ,SINBJ/ 4-1) 

GO  TO  (840,810,820,830),KSYM 
C 

810  XJ»2.*XSYM-XM(J) 

CALL  OC®’F(XP/YP,XJ,YJ,AJ,COSBJ,-SINBJ/-l) 

GO  TO  840 
C 

820  YJ»2.*YSYM-YM(J) 

CALL  OOEFF(XP, YP,XJ, YJ,AJ, -CQSBJ, SINBJ,-1) 

GO  TO  840 
C 

830  XJ«2.*XSYM-XM(J) 

CALL  OOEFF(XP,YP,XJ,YJ,  AJ,COSBJ,~8INBJ,-l) 

XJ-XM(J) 

YJ-2.*YSYM-YM(J) 

CALL  OOEFF(XP, YP,XJ, YJ,AJ,-C06BJ,SINBJ,-1 ) 

XJ*2 . *XSYM-XM( J ) 

CALL  OOEFF(XP, YP,XJ, YJ,  AJ, -CQSBJ, -S1NBJ,+1 ) 

C 

840  CONTINUE 
C 

UX*UXHJXS*D(  JS)+UXN*D(  JN) 

UY-UY+UYS*D  ( JS )  +UYN*D  ( JN ) 

SIGXX>SIGXX+SXXS*D( JS)+SXXN*D( JN) 

SIGYY»SIGYY+SYYS*D<  JS )  +SYYN*D(  JN) 

SI GXY»SIGXY+SXYS*D ( JS ) +SXYN*D ( JN ) 

C 

880  CONTINUE 
C 

W>OINT«NPOINT+l 
IF  (KTYPE.EQ.l)  GO  TO  885 

WRITE  (6,20)  NPOINT,XP,YP,UX,UY,SIGXX,SIGYY,SIGDCY 
GO  TO  890 

885  RAI>9QRT(  (SIGXX-SIGYY)*(SIGXX-SIGYY)/4.+SIGXY*SIQD(Y) 

ANG"PI  -ATAN2  ( -2 .  *8IGXY,  SIGXX-SIGYY )  -2 .  *ANGU 
SIG2«(SIGXX+SIGYY)/2.+RAD*COS(ANG) 
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SIG12*-RAD*SIN(ANG) 

SIFl«SQRT( 2 . *PI *R ) *SlG2 
SIF2=SQRT( 2 . *PI *R ) *SIG12 

WRITE  (6,20)  NPOINT,XP,YP,UX,UY, SIGXX, SIGYY, SIGXY, SIFl,SIF2 
C 

890  CONTINUE 
C 

900  CONTINUE 
C 

C  FORMAT  STATEMENTS. 

C 

1  FORMAT  ( 20A4 ) 

2  FORMAT  <1H1,/,25X,20A4,/) 

3  FORMAT  (314) 

4  FORMAT  (F6.2,E11.4,2F12.4) 

5  FORMAT  (3E11.4) 

6  FORMAT  (/,109H  NUMBER  OF  STRAIGHT-LINE  SEGMENTS  (EACH  CONTAINING  k 
IT  LEAST  ONE  BOUNDARY  ELEMENT)  USED  TO  DEFINE  BOUNDARIES  -,I3,//,12 
23H  NUMBER  OF  STRAIGHT-LINE  SEGMENTS  USED  TO  SPECIFY  OTHER  LOCATION 
3S  (I.E.,  NOT  ON  A  BOUNDARY)  WHERE  RESULTS  ARB  TO  BE  FOUND  -,I3) 

7  FORMAT  (/, 32H  NO  SYMMETRY  CONDITIONS  IMPOSED.) 

8  FORMAT  (/,18H  THE  LINE  X  «  XS  «,F12.4,23H  IS  A  LINE  OF  SYMERRY. ) 

9  FORMAT  (/,18H  THE  LINE  Y  >  YS  «,F12.4,23H  IS  A  LINE  OF  SYK4EIRY. ) 

10  FORMAT  (/,19H  THE  LINES  X  -  XX  -,F12.4,13H  AM)  Y  -  YS  -,F12.4, 

1  >23H  ARE  LINES  OF  SYWETRY. ) 

11  FORMAT  (/ ,18H  POISSON'S  RATIO  -,F6.2,//,18H  YOUNG'S  MODULUS  -,B11. 
14) 

12  FORMAT  (/,31H  XX-OOtPONENT  OF  FI&D  STRESS  -,KU .  4,//,3lH  YY-COWO 
1NBIT  OF  FIELD  STRESS  »,E11 . 4,//, 31H  XY -COMPONENT  OF  FIELD  STRESS  - 
2,  Ell. 4) 

13  FORMAT  (1H1,/,27H  BOUNDARY  EMIT  DATA.,//,96H  ELBCNT  K 

10DE  X  (CENTER)  Y  (CENTER)  LENGTH  ANGLE  UB  OR  SIGHA-S 

2  UN  OR  SIGMA-N,/) 

14  FORMAT  (I4,4F12.4,I4,2E11.4) 

15  FORMAT  ( 219, 3F12 . 4 ,F12 . 2, 2B15 . 4 ) 

16  FORMAT  (1H1,/,66H  DISPLACHMBNTS  AND  STRESSES  AT  MIDPOINTS  OF  B 

10UNDARY  ELEMENTS., //,40H  ELEMENT  D6  US(-)  U8(  +  ), 

2  60H  DN  UN(-)  UN(+)  UX(-)  UY(-)  UX(*), 

3  30H  UY(  +  )  SIGMA-S  SIGMA-N,/) 

17  FORMAT  (I10,10F10.6,2F10.1) 

18  FORMAT  (1H1,/,63H  DI PLACEMENTS  AND  STRESSES  AT  SPECIFIH)  POINT 

IS  IN  THE  BODY.  ,//,117H  POINT  X  CO-CRD  Y  OO-ORD  U 

2X  UY  SIGXX  SIGYY  SIGXY  KI 

3  KII,/) 

19  FORMAT  (4F12. 4,214) 

20  FORMAT  (I9,2F12.4,2F12.6,5F12.1) 

22  FORMAT  (1H1,/,63H  DISPLACEMENTS  AND  STRESSES  AT  SPECIFIED  POINT 
IS  IN  THE  BODY.  ,//,93H  POINT  X  OO-ORD  Y  OO-ORD  UX 

2  UY  SIGXX  SIGYY  SIGXY,/) 

C 

CLOSE  (5) 
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CLOSE  (6) 

END 

SUBROUTINE  INITL 
C 

OOMCN/S2/SXXS/SXXN,SYYS/SYYN, SXYS, SXYN/UXS,UXN,UYS/UYN 
C 

SXXS=0. 

SXXN-O. 

SYYS*0. 

SYYN-0. 

SXYS«0. 

SXYN«0. 

C 

UXS-0. 

UXN-O. 

UYS-0. 

UYN-O. 

C 

RETURN 

END 

SUBROUTINE  CCEEF(X, Y,CX,CY, A,OOSB, SINB,MSYM) 

C 

COUCN/Sl/PI  ,PR/PRl,PR2,OON,OON8 

O0tt10N/S2/SXXS,SXXN,SYYS,SYYN,SXYS,SXYN,UXS,UXN,UYS,UYN 

C 

C0S2B-C06B*C06B-SINB*SINB 
SIN2B-2 . *SINB*OQSB 
C06B2-C06B*C06B 
SINB2»SINB*SINB 

C 

XB*  ( X-CX ) *OOSB+  ( Y-CY  ) *SINB 
YB»-(X-CX)*SINB+( Y-CY)*OOSB 
C 

R1S»(XB-A)*(XB-A)+YB*YB 
R2S- ( XB+A ) * ( XB+A ) +YB*YB 
FL1-0.5*ALOO(R1S) 

PL2»0.5*ALOG(R2S) 

FB2"OON*  ( FL1-FL2 ) 

IF  (YB.NB.O.)  OO  TO  10 
FB3*0. 

IF  (ABS(XB).LT.A)  FB3«C0N*PI 
GO  TO  20 

10  FB3— OCN*(ATAN(  ( XB+A  )/YB )  -ATAN  ( (XB-A)/YB) ) 

20  FB4*OCN*  ( YB/R1S-YH/R2S ) 

FB5-OCN*(  (XB-A)/R1S-(XB+A)/R2S) 

FB6-OON*(  ( (XB-A)**2-YB*YB)/RlS**2-(  (XB+A)**2-YB*YB)/R28**2) 
FB7-2 .  *C0N*YB*  ( ( XB-A )  /Rl8**2-  ( XB+A )  /R2S**2 ) 

C 

UXD6— PR1*8INB*FB2+PR2*C08B*FB3+YB*  ( SINB*FB4-OOSB*FB5 ) 
UXEN— PR1*OOBB*FB2-PR2*SINB*FB3-YB*  (OOSB*FB4+SINB*FB5 ) 
UYD8«PRl*COBB*FB2+FR2*8INB*FB3-YB*  ( C0SB*FB4+8INB*FB5 ) 
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UY DN«- PR1*SINB*FB2+PR2*OOSB*FB3-YB*<SINB*FB4-CQ6B*FB5) 

SXXDS-CCNS*  ( 2 .  *CDSB2*FB4+SIN2B*FB5+YB*  ( OOS2B*FB6-8IN2B»FB7 ) ) 
SXXDNMXNS*  ( -FB5+YB*  ( SIN2B*FB6+C0S2B*FB7 ) ) 

SYYDS-CCNS*  ( 2 .  *SINB2*FB4-SIN2B*FB5-YB*  (COS2B*FB6-SIN2B*FB7 ) ) 
SYYDN-OCNS*  ( -FB5-YB*  ( SIN2B*FB6+CQS2B*FB7 ) ) 

SXYDS-CONS*  ( SIN2B*FB4-CQS2B*FB5+YB*  ( SIN2B*FB6+OOS2B*FB7 ) ) 
SXYDN-CCNS*  ( -YB*  ( COS2B*FB6-SIN2B*FB7 ) ) 

UXS»UXS+MSYM*UXDS 

UXN-UXN+UXDN 

UYS=UYS+MSYM*UYDS 

UYN-UYN+UYDN 

SXXS-SXXS+MSYM*SXXDS 

SXXN-SXXN+SXXDN 

SYY8«8YYS+MSYM*SYYDS 

SYYN-SYYN+SYYDW 

SXYS-SXYS+MSYM*SXYDS 

SXYN-SXYN+SXYDN 

RETURN 

END 

SUBROUTINE  SQLVE(N) 

aXtBN/S3/A<  1600, 1600  ),B(  1600 ) ,  X(  1600 ) 

NB-N-1 

DO  20  J-1,NB 
L-J+l 

DO  20  JJ-L,N 
XM-A(JJ,J)/A(J,J) 

DO  10  I«J,N 

10  A(JJ,I)-A(JJ,I)-A(J,I)*XM 
20  B(JJ)«B(JJ)>B(J)*XM 

X(N)-B(N)/A(N,N) 

DO  40  J«1,NB 
JJ-N-J 
L-JJ+1 
SUM-0. 

DO  30  I-L,N 

30  SUH-8UH+A(JJ,I)*X(I) 

40  X(JJ)«(B(JJ)-SUM)/A(JJ,JJ) 


RETURN 


Appendix  D:  Program  TWODDO 


NOTE:  Most  of  the  Information  in  this  appendix  is  taken  directly 
from  Crouch  and  Starfield  (Reference  7)  with  only  slight  modification. 

1.  Input  (Data)  Deck 

Program  TWODDO  requires  the  following  six  sets  of  cards: 

Set  1:  Free  FORMAT;  one  card  must  be  given  (may  be  blank). 

columns  1-80  of  this  card  contain  any  desired  information  to 
Identify  the  problem  being  solved. 

Set  2:  FORMAT  (314);  one  card  must  be  given  to  specify  the  following 
control  parameters. 

numbs  3  number  of  straight  line  boundary  segments  (each 
containing  at  least  one  boundary  element)  used  to 
define  boundary  contours. 

NUMQS  -  n tuber  of  other  line  segments  (not  art  a  boundary) 
along  which  displacements,  stresses,  and  stress 
.ntensity  factors  are  to  be  computed. 

'  1  no  synmetry  conditions  imposed. 
ksym  ■  2  XkXSYM  (Card  3)  is  a  line  of  symmetry. 

3  y-YSYM  (card  3)  is  a  line  of  symmetry. 

4  x-XSYM  and  y*Y8YM  (card  3)  are  lines  of  symmetry. 

Set  3:  FORMAT  (3E11.4,F6.2,2F12.4);  one  card  must  be  given  to  define 
the  elastic  constants  and  specify  the  locations  of  lines  of 
symmetry  (if  any). 

Ell  *  Young's  modulus  in  x -direction  (Ex ) . 

B22  »  Young's  modulus  in  y-directlon  (By). 

012  ■  Shear  modulus  for  xydlrectlon  (Gxy). 

V12  ■  Poisson's  ratio  for  xy-direction  (vxy). 

XSYM  >  location  of  line  of  symmetry  parallel  to  y-axis  (XSYM 
is  ignored  if  KSYM  *  1  or  3  on  card  2). 

Y8YM  ■  location  of  line  of  symmetry  parallel  to  x-axis  (YSYM 
is  ignored  if  KSYM  *  2  or  4  on  card  2). 

Set  4:  FORMAT  (3E11.4);  one  card  must  be  given  to  define  the  initial 
stresses  (if  any)  in  the  region  of  Interest. 

PXX  ■  oxx  at  infinity. 

PYY  -  o  at  infinity, 
yy 

PXY  ■  r  at  infinity. 

My 


D-l 


Set  5: 


FORMAT  (I4,4F12.4,I4,2E11.4);  NUMBS  cards  must  be  given  to 
define  the  locations  and  boundary  conditions  of  the  boundary 
elements. 


num  ■  number  of  equally  spaced  boundary  elements  along  a 
straight  line  segment,  all  elements  having  the  same 
boundary  conditions. 

XBEG  »  x-coordlnate  of  beginning  of  line  segment. 

YBEX3  -  y -coordinate  of  beginning  of  line  segment. 

XBfl)  ■  x-coordlnate  of  end  of  line  segment. 

YB4D  ■  y -coordinate  of  end  of  line  segment. 


f  1  o  and  of  prescribed. 

2  u*  and  un  prescribed. 

3  u*  and  o'  prescribed. 

4  o'*  and  u^  prescribed. 

BVS  »  resultant  shear  stress  (a  )  or  displacement  (u  ). 

•  • 

BVN  »  resultant  normal  stress  (o's)  or  displacement  (u^). 


Set  6:  format  (4F12.4,2I4);  NUHOS  cards  must  be  given  to  define 
locations  of  points  inside  the  region  of  interest  where 
displacements,  stresses,  and  stress  intensity  factors  are  to  be 
computed. 

xbeg  ■  x-coordlnate  of  first  point  on  line. 

YBEG  »  y -coordinate  of  first  point  on  line. 

XEND  «  x -coordinate  of  last  point  on  line. 

TEND  ■  y -coordinate  of  last  point  on  line. 
nuk>b  >  number  of  equally  spaced  points  between  the  specified 
first  and  last  points. 

KTYPE  >  /  0  compute  displacements  and  stresses  only. 

\  1  compute  stress  Intensity  factors  as  well. 

Note:  if  KTYPE-l,  (XBBG,YEEG)  must  be  the  location  of 
a  crack  tip,  and  the  line  defined  by  (XSMD,yHMD)  must 
extend  In  the  direction  of  the  crack  in  order  far  the 
computed  stress  intensity  factors  to  have  any  meaning. 
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2.  Program  TWODOO  Listing 
PROGRAM  TWODOO 

c 

OOWCW/S1/PI ,  CONI ,  C0N2 ,  G1 ,  G2 ,  ai ,  C22 ,  C12 ,  G12 
OOff10N/S2/SXXS ,  SXXN ,  S  YYS  ,  SYYN ,  8XYS ,  SXYN  ,  UXS ,  UXN ,  UYS  ,  UYN 
OOmON/83/C(  1600, 1600 ),  B<1600 ),  D(  1600 ) 

C 

DIMENSION  XM( 800 ) , YM( 800 ) , A( 800 )  ,OOSfflTT(  800 ) , SINBET(  800 )  ,KOD(  800 ) 
DIMENSION  TITLE(20) 

C 

OPEN  (ONIT-5, FILE* 'TWOTOO. DAT', STATVS-'OLD* ) 

OPEN  (UNIT-6, FILE-' TWODOO.  LI  S',  STATUS-' NEW') 

READ  (5,1)  (TIH£(  I),  1-1,20) 

WRITE  (6,2)  (TITLE( I), 1-1,20) 

READ  (5,3)  NUMBS, NUM06,KSYM 
READ  (5,4)  E11,E22,G12,V12,XSYM,YSYM 
READ  (5,5)  PXX,PYY,PXY 
WRITE  (6,6)  NUMBS, NUM06 
GO  TO  (80,85,90,95),KSYM 
80  WRITE  (6,7) 

GO  TO  100 

85  WRITE  (6,8)  XSYM 
GO  TO  100 

90  WRITE  (6,9)  YSYM 
GO  TO  100 

95  WRITE  (6,10)  XSYM,YSYM 

C 

100  CONTINUE 

WRITE  (6,11)  E11,B22,G12,V12 
WRITE  (6,12)  PXX,PYY,PXY 
C 

PI«4.*ATAN(1.) 

OON-1 .  -V12*V12*E22/E11 
C0N1-E11*G12 

OON2-V12*B22*(  V12*E22/OON+2 .  *G12)-Ell*E22/OCN 
G0N3-B22*O12 

Gl-SQRT(  ( -C0N2+SQRT(  OON2*GON2-4 .  *GCNl*C0N3 )  )/2./O0Nl) 

G2-8QRT( ( -CON2-9QRT(CCN2*OON2-4 .  *OON1*CON3 ) )/2./O0Nl) 

Cll-Ell/OON 

C22-B22/C0N 

C12-V12*C22 

OONl-(Cll*Gl*Gl-G12  )/(C124G12 ) 

OON2-Gl*G2/2 .  /PI/  (C0Nl*G2-Gl/a0Nl ) 

C 

C  DEFINE  LOCATIONS,  SIZES,  ORIENTATIONS  AND  BOUNDARY  CONDITIONS  OF 
C  BOUNDARY  ELB4BNTS. 

C 

NUMBB-0 

DO  110  N-l, NUMBS 

READ  (5,14)  NUM,XBBG,YBBG,XEND,YEND,K00B,BV8,BVN 
XD- ( XEND-XBEG ) /NUM 
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YD-(YEND-YBBG)/»JM 
SW-SQRT  ( XD*XD+YD*YD ) 

C 

DO  110  NE»1,HJM 

NUMBB»NUMBE*1 

M-NUMBB 

XM(M)-XBEG+0.5*(2.*NE-1. )*XD 

YM(M)«YBBG+0.5*(2.*NE-1. )*YD 

A(M)«0.5*SW 

SINBET(M)»YD/SW 

COSBET(M)«XD/SW 

KOD(M)-KODE 

MS«+W-1 
B(MS)*BVS 
110  B(MI)-BVM 
C 

WRITE  (6,13) 

DO  115  M-1,NUMBE 
SIZE*2.*A(M) 

ANOLB>180.*ATAN2(SINBET(M)  ,C0SEET(M)  )/PI 
WRITE  (6,15)  M,KOD(M),  XM(M),  ¥M(M),  SIZE,  ANGUS,  B(2*M-1),B(2*M) 
115  OCNTINUB 
C 

C  ADJUST  STRESS  BOUNDARY  VALUES  TO  ACCOUNT  FOR  INITIAL  STRESSES. 

C  * 

DO  150  N*l, NUMBS 
NN-2*N 
N84W-1 
OOSB-OOSBBT(N) 

SINB-SINBET(N) 

8IG8*(PYY-PXX)*SINB*G0BB+PXY*(G0BB*G08B-8INB*SINB) 
SIGN-PXX*SINB*SINB-2 .  *PXY*SINB*a08B+PYY*008B*C08B 
GO  TO  (120,150,130,140), KOD(N) 

120  B(N8)»B(N8)-SIGB 
B(NN)*B(NN)-SIGN 
GO  TO  150 

130  B(IM)»B(NN)-SIGN 
GO  TO  150 

140  B(N8)*B(NS)'8IGS 
150  CONTINUE 
C 

C  COMPUTE  INFLUENCE  COEFFICIENTS  AND  SET  UP  SYSTEM  OF  ALGEBRAIC 
C  EQUATIONS. 

C 

DO  300  1*1, NUMBS 
IN*2*I 
IS-IN-1 
XI-XM(I) 

YI-YM(I) 

COSBI-COSBETd) 

SINBI-SINBET(I) 
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KODE-KOO(I) 

C 

DO  300  J-1,NUMBE 

JN-2*J 

JS-JN-1 

CALL  INITL 

XJ*XM(J) 

YJ«YM(J) 

COSBJ-COSHBT(J) 

SINBJ-SINBET(J) 

AJ»A(J) 

CALL  OaSFF(XI,YI,XJ,YJ,AJ,COSBJ,SINBJ,+l) 

GO  TO  (240,210,220,230), K8YM 
C 

210  XJ«2.*XSYM-XM(J) 

CALL  OOEFF(XI,YI ,XJ,YJ, AJ,O06BJ,-SINBJ,-l) 

GO  TO  240 
C 

220  YJ«2.*YSYM-YM(J)  . 

CALL  OOHFF(XI,YI,XJ,YJ,AJ,-COSBJ,SINBJ,-1) 

GO  TO  240 
C 

230  XJ*2.*XSYM-XM(J) 

CALL  OOEFF(XI ,YI ,XJ,YJ,AJ,C0SBJ,-SINBJ,-1) 

XJ-XH(J) 

YJ»2.*YSYM-YM(J) 

CALL  OQBFF(XI,YI,XJ,  YJ,AJ,-OQ8BJ,SIMBJ,~l) 

XJ«2.*XSYM-XM(  J) 

CALL  OOHYtXI,  YI,XJ,YJ,AJ,-C06BJ,-8INBJ,+1) 

C 

240  CONTINUE 

GO  TO  ( 250, 260, 270, 280 ) ,KOCG 
C 

250  C(  IS, J8)*( SYY8-SXXS ) *8 I NBI *006BI +8XYS* ( GOBBI *COBBI -SINBI *SINBI ) 
C  ( I S ,  JN )  -  ( SYYN-SXXN ) *SINBI »OOSBI +SXYH* ( OOSBI *C08BI -8IHBI *SIHBI ) 
C(  IN, J8 )«8XXS*SINBI*SINBI-2 . *SXY8*SI NBI *COSBI +SYYS*OQSBI  *OOBBI 
C(  IN,  JN)  -SXXN*8INBI*SINBI  -2 .  *8XYN*SINBI  *C08BI  4-  SYYN*008BI  *C08BI 
GO  TO  300 
C 

260  C(IS,JS)-UXS*G0BBI+UY8*SINBI 
C(  IS,  JN)-UXN*006BI+UYN*SINBI 
C(  IN,  JS)—UXS*8INBI+UYS*C08BI 
C(IN,JN)«-UXN*SINBI+UYN*006BI 
GO  TO  300 
C 

270  C(IS,JS)-UXS*006BI+UY8*SINBI 
C(IS,JN)-UXN*C08BI+UYN*8INBI 

C(  IN, JS ) «8XX8*SINBI *SINBI -2 . *SXYS*SINBI *G08BI+8YYS*C06BI *C06BI 
C( IN , JN ) ■SXXN*SINBI *SINBI -2 . *SXYN*SINBI *00681 +8YYN*008BI *006BI 
GO  TO  300 
C 

280  C ( 1 8, J8 ) - ( 8YY8-8XXS ) *81 NBI »C08BI +8XY8* ( C08BI *C08BI  -8INBI *8INBI ) 


D-3 


C(  I S ,  JN )  -  ( SYYN-SXXN )  *SINBI  *OOSBI +SXYN*  (OOSBI  *CC«BI -8IHBI  *SINBI ) 
C( IN, JS ) — UXS*SINBI+UYS*COSBI 
C( IN, JN )«- UXN*SINBI+UYN*COSBI 
C 

300  CONTINUE 
C 

C  SOLVE  SYSTEM  OF  ALGEBRAIC  EQUATIONS. 

C 

N«2*NUMBE 
CALL  SOLVE(N) 

C 

C  COMPUTE  BOUNDARY  DISPLACEMENTS  AMD  STRESSES. 

C 

WRITE  (6,16) 

DO  600  I«1,NUMBE 
IN«2*I 
IS*IN-1 
XI-XM(I) 

YI»YM(I) 

COSBI-GOSBET(I) 

SINBI-SINBET(I) 

C 

UXNEOO. 

UYNBG-0. 

SIGXX-PXX 

SIGYY-PYY 

SIGXY=PXY 

C 

DO  570  J-1,NUMBE 

JN-2*J 

J8«JN-1 

CALL  INITL 

XJ-XM(J) 

YJ«YM(J) 

AJ-A(J) 

OOSBJ-COSBET(J) 

SINBJ-SINBBT(J) 

CALL  OOffF(XI,YI,XJ,YJ,AJ,COSBJ,SIMBJ,+l) 

00  TO  ( 540, 510, 520, 530 ) ,KSYM 
C 

510  XJ-2.*X8YM-XM(J) 

CALL  GGEFF(XI,YI,XJ,YJ,AJ,OOSBJ,-SINBJ,-1) 

00  TO  540 
C 

520  YJ-2.*YSYM-YM(J) 

CALL  OOEFF(XI,YI,XJ,YJ,AJ,-COSBJ,SINBJ,-1) 

00  TO  540 
C 

530  XJ-2.*XSYM-XM(J) 

CALL  OOEFF (XI,YI,XJ,YJ,AJ, C06BJ , -SI NBJ , -1 ) 

XJ-XM(J) 

YJ-2.*YSYM-YM(J) 
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non 


CALL  C0EFF(XI,YI,XJ,YJ,AJ,-C0SBJ,SINBJ,-1) 

XJ»2.*XSYM-XM(J) 

CALL  COEFF(XI  ,YI,XJ,  YJ,AJ,-G0SBJ,-SINBJ,+1) 

C 

540  CONTINUE 
C 

UXNBG-UXNBG+<JXS*D(  JS ) +UXN*D(  JN ) 

UYNBG-UYNEG+UYS*D(  JS  )HJYN*D(  JN) 

SIGXX-SIGXX+SXXS*D(  JS )  +SXXN*D(  JN ) 

SIGYY-S IGYY+SYYS*D( JS )  +SYYN*D  ( JN ) 

SIGXY-SIGXY+SXYS*D ( JS ) +8XYN*D ( JN ) 

C 

570  CONTINUE 
C 

USNBGHJXNBG»C08BI  +UYNEG*SINBI 
UNNBG--UXNBGPSINBI  HJYNEG*COSBI 
USPOS-USNBG-D(IS) 

UWOS-UNNBG-D(IN) 

UXPOS-USPOS*COSBI -UNPOS*SINBI 
UYPOS-USPOS*SINBI+UNPOS*COSBI 

SIG8-  ( SIGYY-SIGXX )  *SINBI  *CQSBI  +SIGDCY*  ( OQSBI  *CQSBI  -SINBI  *SINBI ) 
SIGN-SIGXX*SINBI*SINBI-2.*SIGXY*SINBI*C0SBH-SIGYY*C0SBI*O06BI 
C 

WRITE  (6,17)  I,D(IS)  ,USNE3G,USPOS,D(IN) ,UMIB8,UW0B,UXMB,UYMBB, 
1  LFXPOS,UYPOS,SIGS,SIGN 

C 

600  CONTINUE 

COMPUTE  DISPLACEMENTS  AND  STRESSES  AT  SPECIFIED  POINTS  IN  BODY. 

IF  (NUM06.LB.0)  GO  TO  900 
NPOINT-O 
DO  900  N-1,NUH0S 

READ  (5,19)  XBBG,YBBG,XEND,YEND,NUMPB,KTYPB 
IF  (KTYPE.BQ.O)  WRITE  (6,22) 

IF  (KTYPE.BQ.l)  WRITE  (6,18) 

NUMP*NUH>B+1 
DELX-(XBND“XBBG)/NWf» 

DELY»(YB(D-YBEG)/NUM? 

IF  (MM>B.OT.O)  NUM>-NUK>+1 
IF  (EGLX**2+D8LY**2.BQ.O. )  MM»1 
C 

ANGLE-ATAN2  ( YEND-YBBG,  XEND-XEBG ) 

C 

DO  890  Hl-l,NUIf> 

XP»XBB0+(NI-1)*DELX 
YP-YBBG+(NI-1 )  *DELY 
C 

UX-O. 

UY-O. 

SIGXX-PXX 

8IGYYM>YY 
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SIGXY-PXY 


C 

00  880  J-1,NUMBE 

JN»2*J 

JS=JN-1 

CALL  INiTL 

XJ*XM(J) 

YJ-YM(J) 

AJ-A(J) 

C 

IF  (SC3RT((XP-XJ)**2+(YP-YJ)**2).LT.2.*AJ)  GO  TO  890 
IF  (KTYPE.BQ.l)  R*SQRT(  (XP-XBBG)*(XP-XBBG)+(YP-YBBG)*(YP-YBBO) ) 
C 

COSBJ*OOSBET( J) 

SINBJ-SINBET(J) 

CALL  CGEFF(XP, YP,XJ, YJ, AJ,006BJ,SINBJ,+1) 

GO  TO  ( 840,810, 820,830 ) ,KSYM 
C 

810  XJ*2.*XSYM-XM( J) 

CALL  OOEFF(XP,YP,XJ, YJ, AJ,C06BJ,-SIMBJ,-1) 

GO  TO  840 
C 

820  YJ«2.*YSYM-YH(J) 

CALL  OC*FF(XP,YP,XJ,YJ,AJ,-C06BJ,SINBJ,-1) 

GO  TO  840 
C 

830  XJ-2.*XSYM-XM(J) 

CALL  CEEFF(XP,YP,XJ,YJ,AJ,a08BJ,-3INBJ,-l) 

XJ-XM(J) 

YJ»2.*YSYM-YM(J) 

CALL  COH*'F(XP,  YP,XJ,  YJ,  AJ,-OQSBJ,SINBJ,-l ) 

XJ»2.*XSYM-XM( J) 

CALL  OOEFF(XP,  YP,XJ,  YJ,  AJ, -C06BJ,-SINBJ,+1) 

C 

840  CONTI  NIIK 
C 

UX-UX+UXS*D(JS)+UXN*D(JN) 

UY-UY+UY8  *D(  JS ) +UYN*D ( JN ) 

SIGXX«SIGXX+SXXS*D(  JS )  +SXXN*D  ( JN ) 
SIGYY-SIGYY+SYY8*D(J8)+SYYN*D(JN) 

8IGXY*8IGXY+8XY8*D( J8)+BXYN*D( JN) 

C 

880  CONTINUE 
C 

NPOINT-tPOINT+1 
IF  (KTYPE.BQ.l)  GO  TO  885 

WRITE  (6,20)  NPOINT,XP,YP,UX,UY,8IGXX,SIGYY,8IGXY 
GO  TO  890 

885  RAD-SQRT< (SIGXX-SIGYY)*<SIGD0(-SIGYY)/4.+SI<Snf*SIGDCY> 

ANG-PI -ATAN2 ( -2 . *8IGXY, 8IGXX-SIGYY ) -2 . *  ANOUK 
8IG2-(8IGXX+8IGYY)/2.+RAD*C08(ANG) 

8IG12— RAD*SIN  ( ANG) 
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SIF1«SQRT( 2.*PI*R)*SIG2 
SIF2-SQRT( 2 . *PI *R ) *3IG12 

WRITS  (6,20)  NPOINT,  XP,  YP,UX,UY,  SIGXX,SIGYY,8IGXY,SIF1,8IF2 
C 

890  CONTINUE 
C 

900  CONTINUE 
C 

C  FORMAT  STA1QCMTS. 

C 

1  FORMAT  (20A4) 

2  FORMAT  (1H1,/,25X,20A4,/) 

3  FORMAT  (314) 

4  FORMAT  (3E11.4,F6.2,2F12.4) 

5  FORMAT  (3B11.4) 

6  FORMAT  (/,109H  NUMBER  OF  STRAIGHT-LINE  SEGMENTS  (EACH  CONTAINING  A 
IT  LEAST  ONE  BOUNDARY  BJ3MENT)  USED  TO  DEFINE  BOUNDARIES  «,I3,//,12 
23H  NUMBER  OF  8TRAIGHT-LINE  SEGMENTS  USED  TO  SPECIFY  OTHER  LOCATION 
3S  (I.B.,  NOT  ON  A  BOUNDARY)  WWRE  RESULTS  ARE  TO  BE  FOUND  «,I3) 

7  FORMAT  (/,32H  NO  8Y)#f3TRY  CONDITIONS  IMPOSED.) 

8  FORMAT  (/,18H  THE  LINE  X  -  XS  «,F12.4,23H  IS  A  LIMB  OF  SYMMETRY.) 

9  FORMAT  (/,18H  THE  LINE  Y  ■  YS  -,F12.4,23H  IS  A  LINE  OF  SYMMETRY. ) 

10  FORMAT  (/,19H  THE  LINES  X  -  XS  -,F12.4,13H  AND  Y  -  YS  -,F12.4, 

1  23H  ARE  LINES  OF  SYMMETRY. ) 

11  FORMAT  (/, 14H  MODULUS  EXX  «,K1I.4,/,14H  MODULUS  EYY  «,B11.4,/,14H 
‘1MOOULU8  GXY  -,E11.4,//,20H  POISSON  RATIO  VXY  «,F6.2) 

12  FORMAT  (/,31H  XX-OCM’GNBNT  OF  FIBJ>  STRESS  «,E11.4,//,31H  YY-OOWO 
INERT  OF  FIELD  STRESS  «,Kll.  4,//,31H  XY-COWOMNT  OF  FIBJ>  STRESS  * 
2, Ell. 4) 

13  FORMAT  (1H1,/,27H  BOUNDARY  ELEMENT  DATA.,//,96H  HJ94BNT  '  K 

IOCS  X  (CB4TBR)  Y  (CSfRR)  LENGTH  ANGLE  US  OR  SIGMA-S 

2  UN  OR  8IGMA-N,/) 

14  FORMAT  (I4,4F12.4,I4,2R11.4) 

15  FORMAT  (2I9,3F12.4,F12.2,2E15.4) 

16  FORMAT  (1H1,/,66H  DISPLACEMENTS  AND  8TRE8BBB  AT  MIDPOINTS  OF  B 

10UNDARY  ELEMB4TS . , //, 40H  HMNT  DS  US(-)  U8(+), 

2  60H  DN  UN(-)  UN( ♦ )  UX(-)  UY(-)  UX<4-), 

3  3 OH  UY(+)  SIGMA-S  SIOMA-N,/) 

17  FORMAT  (I10,10F10.6,2F10.1) 

18  FORMAT  (1H1,/,63H  DISPLACEMENTS  AND  STRESSES  AT  SPECIFIED  POINT 

IS  IN  THE  BODY . ,  ,  117H  POINT  X  CO-ORD  Y  CO-ORD  U 

2X  UY  8IGXX  8IGYY  SIGXY  KI 

3  KII,/) 

19  FORMAT  (4F12. 4,214) 

20  FORMAT  (I9,2F12.4,2F12.6,5F12.1) 

22  FORMAT  (1H1,/,63H  DISPLACBENTS  AM)  STRESSES  AT  SPECIFIED  POINT 
18  IN  TW  BODY. ,//,93H  POINT  X  CO-ORD  Y  CO-ORD  UX 

2  UY  8IGXX  8IOYY  8IGDCY,/) 

C 

CLOSE  (5) 

CLOSE  (6) 

BID 


SUBROUTINE  INITL 
C 

C0MM0N/S2/SXXS , SXXN , SYYS , SYYN, 8XYS , 8XYN, UX8 , UXN , UYS , UYN 
C 

SXXS-O. 

SXXN»0. 

SYYS-0. 

SYYN»0. 

SXYS*0. 

SXYN-0. 

C 

UXS-O. 

UXN-0. 

UYS-0. 

UYN*0. 

C 

RETURN 

EM) 

SUBROUTINE  OCEFF(X,  Y,CX,CY,  A,COGB,SINB,MSYM) 

C 

OOMCN/S1/PI ,C0Nl,C0N2,Gl,G2,Cll,C22,Cl2,C66 
COMON/S2/SXXS,  SXXN,  SYYS,  SYYN,  SXYS,SXYN,UXS, UXN, UYS,UYN 
C 

CDS2B-CQSB*OOSB-SINB*SINB 
SIN2B-2 .  *SINB*OOSB 
006B2«006B*006B 
SINB2«SINB*SINB 

C 

XB-(X-CX)*O06B+(Y-C5f)*8INB 
YB—  ( X-CX ) *SINB+  ( Y-CY  > *COBB 
C 

Al*  ( Gl*Gl*006B2+8INB2 )  /Gl/Gl 
A2«(G2*G2*C08B2+SINB2  )/G2/G2 
Bl«  ( 1 .  -Gl*Gl )  *SIN2B/G1/G1 
B2«(l.-G2*02)*8IN2B/02/02 
Cl»(Gl*Gl*8INB24O06B2 ) /Gl/Gl 
C2«(G2*G2*SINB2+O06B2  )/G2/G2 
Rl81“Al*(XB-A)*(XB-A)+Bl*(XB-A)*YB*d*YB*YB 
R182-A2* ( XB-A ) * ( XB-A ) +B2* ( XB-A) *YB4C2*YB*YB 
R281«Al*<XBfA)*(XB+A)4>Bl*(XB+A)*YB4d*YB*YB 
R2S2-A2*(XB+A)*(XB+A)+B2*(XB+A)*YB+C2*YB*YB 
IF  (YB.NB.O.)  GO  TO  10 
FB1*0. 

FB2*0. 

IP  (ABS(XB).GT.A)  GO  TO  20 
FBI— PI 
FB2— PI 
GO  TO  20 

10  FB1-ATAN2 ( YB/G1/A1 , XB-A+Bl*YB/Al/2 . ) -ATAN2 ( YB/Gl/Al , XB+A+H1*YB/Al/ 

12.) 

FB2-ATAN2  ( YB/G2/A2,  XB-A+B2*YB/A2/2 . ) -ATAN2  ( YB/G2/A2 ,  XB+A+B2*YB/A2/ 

12.) 
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20  FB3-(AL0G(R1S1/R2S1) )/2. 

FB4-(ALOG(R182/R2S2))/2. 

R11—YB/OI/R1S1+YB/G1/R2S1 
R12— YB/G2/R1S2+YH/G2/R2S2 
R21* ( XB-A)/G1/R1S1-(XB+A)/G1/R2S1 
R22»(XB-A)/G2/RIS2-(XB+A)/G2/R2S2 

R31* ( Al* { XB-A )  +Bl*YB/2 .  )/R181-(Al*(XB+A)+Bl*YB/2.  )/R2Sl 
R32* ( A2* ( XB-A )  +B2*YB/2 .  )/RlS2-(A2*(XB+A)+B2*YB/2.  )/R2S2 
R41»(Bl*(XB-A)/2.+Cl*YB)/R181-(Bl*(XB+A)/2.+Cl*YB)/R281 
R42«  ( B2*  ( XB-A ) /2 . 4C2*YB ) /R1S2- ( B2*  ( XB+A  )/2 . 4C2*YB  )/R2S2 
FB5-R11/OOH1/G2-R12*OOH1/G1 
FB6«R21/OOM1/G2-R22*OCN1/G1 
FB7*R11*OON1/G1-R12/OON1/G2 
FB8*R21*CCN1/Gl-R22/<XN1/G2 

UXDS“+OOH2*  (CX)SB*  (FB1/OOW1/G2-PB2*OON1/G1 )  +8INB*  (FB3-FB4 ) ) 

UXDN— CXW2*(SINB*(FBl/C30HVG2-PB2*O0Hl/Gl)-aMB*(FB3-FB4) ) 
UYDS*+CCH2*  ( OQSB/G1/G2*  ( FB3-FB4 )  -SINB*  ( FB1*OOH1/G1-FB2/OOH1/G2 ) ) 
UYDN-- CXJN2*  ( SINB/G1/G2*  ( FB3-PB4 )  +OQ8B*  ( FB1*OOW1/Gl-FB2/OOH1/G2 ) ) 

SXXDS-+OON2*  ( dl*COSB*  ( G06B*FB5-SINB*FB6 )  -C12*SIHB*  ( SIKB*FB7+008B* 
1FB8 )  +  (C11+C12/G1/G2 )  *SINB*OOBB*  (R31-R32 )  -(C11*SINB2-C12*OOBB2/G1/G 
22)*(R41-R42) ) 

SXXEN«-CCN2*  (CU*SINB*  ( 008B*FB5-8INB*FB6 )  +0.2*0088*  (  SINB*FB7+OOSB* 
1FB8 )  -  (Cai*OOSB2-C12*SINB2/Gl/G2 )  *  (R31-R32 )  ♦  (Cll+d2/01/G2 )  *SINB*CO 
2SB*(R41-R42) ) 

SYYDS-+CCH2*  (02*0088*  ( G08B*FB5-8INB*FB6  )  -C22*SINB*  ( SINB*FB7+COSB* 
IFB8 )  +■  ( Cl  2  +C2  2/G1/G2 )  *SINB*OOBB*  ( R31-R32 )  -  ( d2*8INB2-C22*008B2/GL/G 
22)*(R41-R42) ) 

SYYDH— CXW2*  (CL2*8I  MB*  ( G08B*FB5-SINB*FB6 )  +C22*OOSB*  ( SINB*FB7+O08B* 
1FB8 ) -(O2*C0SB2-C22*SINB2/Gl/G2 )  *  (R31-R32 )  *■  (024022/01/02 )  *8INB*G0 
2SB*(R41-R42) ) 

SXYD8«+OON2*C66*  ( SI NB*00SB*  ( FB5-FB7 )  +C08B2*FB6+SINB2*FB8+ (  SINB2+0D 
1SB2/G1/G2  )*(R31-R32 )  ♦8IHB*C08B*  ( 1 .  -1  ./Gl/02  )*  (R41-R42 ) ) 
SXYEN«-OON2*C66*  ( SINB2*FB54006B2*FB7+S1NB*006B*  (FB6-FB8  )  -8INB*C08B 
1*  ( 1 .  -1./G1/G2 )  *(R31-R32 ) -(O08B2+8INB2/OI/G2  )*(R41-R42 ) ) 

UXS-UX8++I8YM*UXD6 

UXN-UXN-HJXDN 

UY8-UY8+M8YM*UYD6 

UW-UXN4UXDH 

SXX8-8XX84H8YM*8XXD6 

8XXN-8XXN+SXXEN 

SYYS«8YY8+MBYM*8YYD8 

SYYN-SYYN4-8YYEN 

SXYS-SXY84KBYM*8XYD6 

SXYN-8XYN+8XYEN 

RETURN 

END 

SUBROUTINE  SOLVE(N) 


D-ll 


c 


COttW/83/A  ( 1600 , 1600  ),B(  1600 ) ,  X  ( 1600 ) 


NB-N-1 

DO  20  J-1,NB 

L-J+l 

DO  20  JJ«L,N 
XM*A( JJ, J )/A( J, J) 

DO  10  I-J,N 

10  A(JJ,I)«A<JJ,I)-A(J,I)*XM 
20  B(JJ)«B(JJ)-B(J)*XM 

X(N)«B(N)/A(N,N) 

DO  40  J*1,NB 
JJ-N-J 
L-JJ+1 
SUM*0. 

DO  30  I-L,N 

30  SUM-SUM+A(JJ,I)*X(I) 

40  X(JJ)«(B(JJ)-8UH)/A(JJ,JJ) 

RETURN 

END 
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